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1. INTRODUCTION 


Work on the prediction of unsteady airloads on helicopter rotors has been underway 
for several decades, as documented by References 1-8. One motivation for this effort has 
been the high level of vibration experienced by rotorcraft in a variety of flight conditions, 
much of which is caused by unsteady loading on helicopter main rotors. On the whole, it 
can be said the understanding of the dynamic response of rotor blades to aerodynamic 
loading is well advanced, while the prediction of the wake-induced velocities that produce 
this loading is more uncertain. The prediction of the velocity field induced by the rotor 
wake and the wake-induced aerodynamic loads are the primary topics of interest in the 
effort described here. 

It is axiomatic that the study of helicopter rotor aerodynamics in forward flight is an 
exceptionally complex undertaking involving most of the major phenomena being studied 
by modem aerodynamicists, e.g., compressible subsonic and transonic flow, viscous drag, 
stalled and separated flows, and unsteady aerodynamics. Features of all of these 
phenomena should be included in a comprehensive model of rotor aerodynamics. 
However, it was not judged desirable to attempt to include the state-of-the-art treatments of 
all of these issues within this analysis, since such treatments often involve direct integration 
of the primitive equations of fluid mechanics (e.g., the Euler or Navier-Stokes equations) 
and so typically require exceptionally large computational capabilities. Such an approach 
would have defeated the primary purpose of this effort: the development of an analysis of 
rotors in forward flight that not only represents a substantial improvement in the refinement 
of rotor wake modeling but is also usable in the near term as a research tool in support of 
design efforts in the rotorcraft industry. 

As will be made clear in the discussion to follow, the primary technical focus of the 
present work has been the implementation of an advanced mcxlel of the vortex wake for use 
in a general analysis of helicopter rotor aerodynamic loading in forward flight. This model 
features a novel discretization of the wake into filamentary vortices representing constant- 
strength contours of vorticity in the vortex sheets that trail from each blade. This approach 
automatically captures the wake of the full span of the rotor blade and so contains structures 
absent from earlier, more simplified models. The new full-span wake model is also 
noteworthy in that it eliminates the artificial distinction between vorticity generated by 
azimuthal and spanwise changes in the bound circulation on the blade (i.e., between what 
is conventionally termed "trailed" and "shed" vorticity) and in that it provides a visually 
meaningful representation of the wake. 

The wake model itself has been constructed in a sufficiently general fashion to 
capture the full semi-infinite wake of the helicopter rotor, and includes efficient and 
accurate calculation of far wake effects. The computation of aerodynamic loads in the 
presence of the wake-induced flow field is carried out using a vortex lattice model of the 
blade. This model is appropriate for rotor calculations because of its potential for refined 
treatment of tip effects and wake/blade interaction. In addition, the study of helicopter 
blade loading is a classic example of a truly aeroelastic calculation, and so an appropriate 
treatment of the dynamics of the rotor blade is particularly important. It was not the 
intention of the present effort to break new ground in the analysis of rotor blade dynamics, 
so the model developed here includes the implementation of a finite element analysis of 
blade structure drawn largely from the existing literature on structural dynamics. The 
resulting coupled analysis of vortex wake dynamics, unsteady blade loading, and blade 
deformation is suitable for application to isolated rotors in steady forward flight. 



Given the focus of the present work, the Computation of Rotor Aerodynamics in 
Forward Flight, the complete analysis has been dubbed RotorCRAFT. This report 
discusses the fundamental work that was carried out to support the development of the 
initial version of RotorCRAFT and documents the model as it currently exists. It opens 
with a brief review of background information and literature on the prediction of forward 
flight airloads, providing the motivation for the development of the full-span free wake 
model that forms the heart of the RotorCRAFT analysis. The following section discusses 
the principles of the wake model and the major features of its implementation. Along with 
this is a description of the implementation of a new class of vortex wake elements based on 
Analytical Numerical Matching (ANM), elements that offer the promise of substantially 
increased efficiency with no sacrifice in accuracy relative to the Basic Curved Vortex 
Elements (BCVEs) currently implemented in the model. 

The subsequent discussion outlines the application of vortex lattice methods to 
aerodynamic load calculations and the implementation of the finite element model of blade 
deformation. The remainder of the report focuses on describing the structure and operation 
of RotorCRAFT and its performance in a variety of correlation studies which were aimed at 
evaluating the code's success as a computational tool for the analysis of aerodynamic loads. 
The discussion of these results will point out the advantages gained through the application 
of full-span wake modeling to rotor loads prediction and will note the features of the 
analysis requiring further development. 
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2. BACKGROUND ON AIRLOAD PREDICTIONS IN FORWARD FLIGHT 


Miller (Ref. 1) discussed many of the sources and consequences of vibratory 
airloads in helicopters, noting that the higher harmonic components of downwash persist 
even in high speed flight, contrary to the trends predicted by models based on uniform 
downwash. He also pointed out the importance in capturing wake distortion so as to 
properly account for unsteady loading in both low and high speed flight. Scully (Ref. 2) 
developed a free wake analysis of the rotor using two free filaments (a rolled-up tip vortex 
and a diffuse inboard filament) to model the rotor wake. This distorted wake model 
succeeded in capturing certain features of the unsteady loading on rotors in a variety of 
flight conditions, though uncertainties about the proper modelling of close blade/wake 
interactions precluded accurate prediction of higher harmonic loads. Other broadly similar 
wake models have been developed over the last fifteen years (e.g., the work of Egolf and 
Landgrebe, Ref. 3, Sadler, Ref. 4, and Johnson, Ref. 5). 

The seminal paper of Hooper (Ref. 6) pointed out the inadequacy of these types of 
models in predicting the nearly impulsive airloading events that occur in a broad range of 
rotor designs, particularly in high speed flight. Hooper considered data for flight 
conditions at both low and high advance ratios. He concluded that the major features of 
vibratory loads at low forward speed could be predicted using relatively simple models 
involving the oblique interaction of blades with straight line vortices. However, the high- 
speed cases proved to contain features requiring substantially more advanced modeling 
techniques. Hooper showed that several ostensibly dissimilar rotors shared a strikingly 
similar loading pattern and that in several cases the higher harmonic airloading was 
dominated by a relatively discrete event near the blade tip on the advancing side. During 
such interactions the blade experiences a sharp upward, and then downward, load 
fluctuation. This up-down pulse, which is clearly an important contributor to the overall 
higher harmonic airloading, appeared to be due primarily to aerodynamic interaction with 
the wakes produced by previous blades. Hooper also documented the inability of 
conventional rotor wake models to compute the behavior observed in the experimental data, 
a result suggesting that some aspects of the traditional approach to rotor wake 
aerodynamics must be inadequate, at least for the high speed case. Recent work by Miller 
(Ref. 7) and Johnson (Ref. 8) has been directed toward repairing some of the shortcomings 
in previous single-tip-vortex models. 

Any w ake model designed to address the issue of higher harmonic airloading must 
be sufficiently detailed to capture the important details of the loading experienced by rotor 
blades. Figure 2-1 shows that rapid, complicated load variations occur on the advancing 

side. Around vy = 90* there is a region of negative load at the tip, and a relatively gradual 
span wise load variation, with the maximum load shifted far inboard. The advancing blade 
experiences a rapid decrease in loading as it enters this region and a rapid load increase 
upon leaving it. These load variations generate a complex wake field which is encountered 
by the following blades. Strong trailed (streamwise) and shed (spanwise) vorticity 
components are present in this incident wake field. Vorticity of opposite sign is generated 
by the negative loading near the tip of the generating blade, resulting in a negative sign tip 
vortex. The usual positive sign tip vorticity is also present but is shifted inboard. Because 
of the more gradual radial load variation, there is insufficient time for this positive sign 
sheet to co mp letely roll up before encountering the following blades. 

Clearly, in view of this loading pattern and its likely source, a new treatment of the 
wake is necessary if rotor unsteady airloads are to be adequately predicted. The technical 
approach to be described below utilizes a refined analysis of the near wake that extends 
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Figure 2-1. Linear surface plots of aerodynamic load distribution on the 
H-34 main rotor at advance ratio 0.39 (from Ref. 6). 
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beyond the first few interactions with following blades. This analysis predicts the 
convective evolution of the actual wake structure, and the approach is based on calculating, 
rather than modeling, the actual wake flow field. For this reason, the analysis is designed 
to acco mm odate a rather general description of the wake. The wake structure accounts for 
the effects of spatial and temporal load variations on the generating blade. Preliminary 
development of these concepts is described in References 9-17, which contain very 
promising initial results on the predictions of the vibratory airloading of rotors, as well as 
discussion of fundamental tools for rotor wake analysis. The work described in 
References 12 and 13, which originated as part of the Boeing Helicopters research 
program, was of particular importance in demonstrating the applicability of the general 
approach used here to the prediction of rotor airloads in high speed flight. 

It is important to note, though, that despite its importance to the prediction of 
helicopter vibration, the issue of high frequency airloading cannot be considered in 
isolation. The complete aerodynamic load calculation on the rotor - including the coupling 
to performance, trim, and blade dynamics calculations - must be considered as a whole. 
For this reason, this report goes considerably beyond previous exploratory work and 
addresses integrated and low-frequency rotor loading in substantial detail, as is appropriate 
for a general analysis of rotor aerodynamics. Indeed, though Hooper's discussion in 
Reference 6 correctly notes that the fundamental mechanisms leading to vibratory loads in 
transition flight at low advance ratio have been reliably identified, the predictive capability 
of current rotor aerodynamic models in this regime is at best fair, as was made clear by a 
comprehensive review paper by Harris (Ref. 18) as well as more focused studies of low- 
speed rotor behavior (e.g.. Ref. 11). Thus, the discussion that follows will emphasize the 
full spectrum of rotor loads calculations and will include correlation studies with rotors in a 
broad range of flight conditions, in addition to the high-speed cases described above. 
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3. FULL-SPAN FREE WAKE MODELING 


This section describes two of the principal computational tools used in the wake 
simulations featured in this report: curved vortex elements and the novel full-span wake 
model. The application of curved elements is briefly reviewed, followed by a discussion of 
the representation used for the complex wakes characteristic of forward flight. The model 
described here is based on the same fundamental ideas as the forward flight wake model 
outlined in Reference 12, but has been radically altered and extended by several major 
subsequent modifications. These modifications, documented in part by References 14 and 
15, include the addition of a vortex lattice blade model, a more consistent wake evolution 
scheme, implementation of a far wake model, and a coupling to a flexible finite element 
structural model, as well as a variety of additional features developed during the current 
effort. 


In addition, the current wake model takes advantage of recent work by Bliss and 
Miller (Refs. 16 and 17) in the development of exceptionally efficient techniques for the 
prediction of vortex-induced flow fields in free wake calculations. This work, which 
involves the application of a technique known as Analytical/Numerical Matching (ANM), 
has resulted in the development of methods for the prediction of vortex-induced velocity 
fields that are in many ways superior to the curved element model implemented in the 
earlier versions of the RotorCRAFT analysis. While not yet sufficiently mature to replace 
wake models based on curved elements, the ANM wake model is available as an option 
within RotorCRAFT. For completeness, the fundamental features of the entire free wake 
flow field model are now briefly reviewed, including discussion of curved vortex elements, 
full-span wake modeling, and ANM. 


3.1 Curved Vortex Elements in Free Wake Models 

The purpose of a free wake analysis is to simulate in detail the actual shape and 
motion of the wake of each rotor blade treated as a free vortex flow. The wake is 
represented by vortex filaments trailed from the blades, and each trailing filament is 
composed of a string of individual vortex elements connected at collocation points. The 
velocity field of these individual vortex elements is determined through numerical 
evaluation of analytical expressions for the velocity field of these elements. The Biot- 
Savart integration for the overall wake velocity field is evaluated numerically by summing 
the contributions of the individual elements. In conventional free wake analyses like that 
described in References 10 and 11, the convection of these points over a small time 
interval, corresponding to an increment of rotor rotation, is determined at each step given 
this velocity field. The vortex elements are then repositioned between the new collocation 
point locations and the calculation procedure is repeated. The process can be considered 
converged when the results are repeatable from one rotor revolution to the next 

Two types of vortex elements are required to implement a free wake analysis. The 
first type of element is used to evaluate the velocity induced at any point in the flow field, 
except points on the element itself. Traditionally, straight-line elements have been used for 
this purpose (see for example Refs. 2-5). An alternate approach was developed using a 
more sophisticated curved element (Refs. 9-11). This element is referred to as the Basic 
Curved Vortex Element, or BCVE. A second type of element is used only to evaluate the 
velocity induced by the element on itself. This element is called the Self-Induction Vortex 
Element, or SIVE. Curved elements have previously been applied to work on the wakes of 
hovering rotors as well as simulations in low- and moderate-speed forward flight (Ref. 
11). Reference 10 also contains a discussion of the gains in efficiency and accuracy that 
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are realized when curved elements are used in place of straight-line elements. The 
documented success of curved elements in wake modeling is an important facet of the 
forward flight wake model described in this report. 

Figure 3-1 shows a typical BCVE, depicting the local coordinate system associated 
with it as well as its parabolic shape. For such a parabolic shape, the functional form of the 
Biot-Savart law that is integrated over its length to compute the velocity field is such that it 
can be solved analytically; this circumstance permits a very favorable tradeoff to be made in 
accuracy and computation time when calculations using BCVEs are compared to ones 
involving straight line elements. Figure 3-2a shows several methods for joining BCVEs 
together to form the continuous vortex filaments that are the building blocks of most current 
Lagrangian free wake models. At present, the "interpolated point" method is used, and the 
geometry of the filaments at any instant is determined by passing overlapping arcs between 
the existing collocation points. Also evident in Figure 3- 2b is the circular arc SIVE that 
computes die velocity field in the immediate vicinity of the collocation points themselves. 
The construction of vortex filaments using BCVEs and SIVEs is discussed in detail in 
Reference 10. 


3.2 Representation of Full-Span Wakes 

This section outlines the full-span wake model that forms the heart of the 
RotorCRAFT analysis. As discussed in Section 2, a complex wake model is needed in 
forward flight since the time- varying loading experienced by rotor blades generates both 
trailed and shed vorticity components in the wake. For reference, note the typical 
experimental result for aerodynamic load variations that are experienced by a rotor blade in 
high speed flight as shown in Figure 2-1. An idealized sequence of advancing side load 
distributions containing this behavior is sketched in Figure 3-3. A corresponding sequence 
of bound circulation distributions would exhibit similar behavior. Spanwise variations in 
the bound circulation are responsible for the "trailed" component of wake vorticity and 
temporal (Le, azimuthal) variations are responsible for the "shed" component. 

Assume that the bound circulation distribution along the span of the rotor at any 

azimuth angle \|/ can be represented as r(r,y), where r is the radial distance along the blade 
span. The wake that trails from the blade consists of a continuous sheet of vorticity whose 

strength as it leaves the blade is y(r,\j/). The vorticity in the wake is a vector quantity, and 
can be represented by 

Y = Yx i + YyJ (3.1) 

where x and y are defined parallel to the span and to the chord of the blade, respectively. 
The intensity of each component of vorticity is related to the bound circulation by 

_ -ar(r,y) 

Yx dt (3.2) 

and 

-ar(r,y) 

Yy at (3.3) 
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END-TO-END CONNECTION METHOD 



OVERLAPPING ARCS METHOD 



OVERLAPPING ARCS 


INTERPOLATED POINT METHOD 


Figure 3-2a. Connection methods for curved vortex elements. 



Figure 3-2b. Typical arrangement of elements to calculate the self-induced 
velocity at point j . 
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Figure 3-3. Typical load variations on the advancing side in 
high speed flight. 
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while the magnitude of the vector is determined by 


Y = f3 = Vy x 2 + Y y 2 (3.4) 

Note that the strength of the wake sheet may be zero and that each component may reverse 
sign depending on the strength and rate of change (spatial or temporal) of the bound 
circulation. 

Figure 3-4 shows the idealized wake vorticity field corresponding to the load 
variations in Figure 3-3. The wake here is shown in terms of contours of constant sheet 

strength, i.e. lines of constant y are depicted. (Note that for the current discussion, the 
distortion of sheet due to wake on wake interaction is ignored and the wake is assumed to 
lie in a plane; this situation is well approximated by a rotor in high speed forward flight). 
Clearly, only a finite number of contour lines are used in this figure, though a continuum of 

lines are in fact required to completely represent the wake sheet. Because increments in y 
are constant between each contour line, the amount of circulation contained between any 
two contour lines on the sheet is constant. This circulation is related to the bound 
circulation on the blade as follows: 


Ar = 



+ Ar 


7 (r,y) dr = r(r+Ar,y) - r(r,y) 


(3.5) 


The contour lines may come together to form a closed loop, but other than this the total 
circulation between any two lines will remain constant throughout the wake. 

One consequence of this is that if the origination points of contour lines (denoted 
"release points” in the discussion below) are spaced radially along the blade corresponding 
to fixed increments in bound circulation, then the radial distance between the contour lines 
will be a direct measure of the gradient in bound circulation along the blade. That is, since 

AT is a constant, the approximate relationship 


ar(r,y) s A r 

^ Ar(r,y) (3.6) 

means that filament spacing will be inversely proportional to span wise circulation gradient. 
(Here, Ar is the radial spacing between contour lines as they leave the blade.) Thus, 
contour lines will be tightly spaced in the radial direction where gradients of bound 
circulation are high, as is typically true near the blade tip at most azimuth angles. As Figure 
3-3 makes clear, however, realistic blade load variations do not always produce large 
bound circulation gradients near the blade tips. Gradients are often shallow at the 
advancing tip in high speed flight, and so the wake may take on a more sheetlike character, 
rather than rolling up into a concentrated tip vortex. Conversely, bound circulation 
gradients may not always be shallow in certain regions. Note in particular that the radial 
gradients near the blade root may be relatively large compared to the weak gradients 
typically seen in the inboard loading distribution of hovering rotors. Thus, unlike the 
diffuse inboard wakes usually present on rotors in hover, concentrated vortical structures 
may be present trailing from the root region on rotors in high speed flight. 
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Note that the tangent to the contour line at any point determines the direction of the 
*♦ 

local vorticity vector 1 . Often, the terms "trailed" and "shed" vorticity are used in 

describing the wake, the former to denote vorticity trailed due to spatial (radial) gradients in 

bound circulation and the latter to denote vorticity due to temporal (azimuthal) gradients. 

As should be clear from the discussion to this point, these two quantities are in fact 

-* 

components of the resultant vector 7 . Clearly, it is desirable to employ a wake 
discretization that removes the artificial distinction between these two quantities and treats 
the vector vorticity field in a unified manner. 

In principle, one could achieve this end by discretizing the wake into continuous 
sheet elements, but this has proved awkward in the past (Ref. 2); curved sheet elements 
that could smoothly capture the distortion of the wake sheet are computationally inefficient, 
while simpler flat sheet elements tend to develop kinks and gaps. Another possible 
representation of the wake sheet is a lattice of discrete shed and trailed vortex elements (see, 
for example, Ref. 4) as shown in Figure 3-5. This lattice is typically constructed by 
trailing mutually perpendicular straight line vortex elements from the blade, elements 
parallel to the blade span to account for the x component of vorticity, while elements 
perpendicular to the span capture the y component. The trailing components have an 

azimuthal extent of A\jt, while the length of the shed components is determined by the radial 
spacing of the vortex lattices on the blade. This approach has been used previously to 
calculate the near wake downwash on the generating blade (Refs. 4, 19, and 20). By 
extending the lattice far enough downstream, the wake interaction with subsequent blades 
can also be simulated. In fact, a lattice representation using curved vortex elements could 
be considered for use in a free wake treatment. However, in a free wake analysis this 
approach seems unnatural, in part because it preserves the artificial distinction between 
shed and trailed vorticity noted above. In addition, the effect of discretized shed 
components on the roll-up of the trailed components is open to question. 

A wake model built on the contours of constant wake vorticity is preferable. 
Several interesting features of the wake vorticity field are apparent, in addition to those 
already cited. First, note that the radial location of the release points of the constant contour 
lines changes in accordance with the varying bound circulation distribution. Also, consider 
the path of the "zero contour" line that trails from the locations of the maximum values of 
bound circulation at each azimuth (the locus of these maxima are shown with a dashed line 
in Figure 3-4). When the magnitude of this maximum value changes by an amount larger 

than the increment of circulation AT between the contour lines, the lines on either side of 
this maximum line must come together to form loops. The loops across the maximum line 
in Figure 3-4 are a consequence of decrease and subsequent increase in the maximum load 
(and bound circulation) sketched in Figure 3-3. 

As well as for min g on either side of the line of local maxima on the rotor blade, 
loops may form around local mi nima as well. Consider the loading pattern sketched in 
Figure 3-3. There, the negative bound circulation that appears and then disappears near the 
tip leads to the formation of the set of elongated closed contours seen in the tip region in 
Figure 3-4. The number of loops is roughly 
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Figure 3-5. Modeling the trailed and sheet vorticity by a lattice 
of vortex elements. 
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The precise integral number of loops trailing from this tip region will of course be a 
function of the discreteness of the contour representation. 

The discussion to this point has focused on the qualitative behavior of vortex sheet 
strength contours trailing from the blade, but has not explicitly defined a method for 
carrying out practical wake calculations. The first step in motivating this model is to recall 
that the contour lines can in fact be treated as vortex filaments. The natural step to take in 
representing the wake is thus to lay down curved vortex elements along the contours of 
constant wake strength, as shown in Figure 3-6. To conserve circulation in the wake, 
these filaments are assumed to contain all of the rolled-up vorticity trailed between adjacent 

release points. This means that all of the filaments will have equal strength AT, and the 
filaments have constant strength along their length. 

The shed and trailed vorticity is thus accounted for by the fact that the direction of 
the tangent vector to the resultant elements changes continuously even though the strength 
of each filament is constant . Using this approach, on the average only half as many 
elements are required as in the straight-line lattice method, where the local vorticity vector is 
represented by two discrete straight elements. Thus, the use of curved elements means that 
fewer, relatively large elements can be used to accurately represent the contour shapes. 
These factors give the full-span free wake model based on curved vortex elements a 
substantial advantage in efficiency relative to the lattice approach; since the model requires 
roughly half the number of elements, the number of operations required to compute the 
wake-on-wake velocity field will be reduced by approximately a factor of four. Previous 
work described in Reference 10 has indicated that the ratio of execution time for wake 
models using curved or straight vortex elements is 



The factor of 3.3 used here was determined by numerical experiments to be the ratio of 
execution time of BCVEs to straight elements (Ref. 10). Clearly, since a single curved 
element can replace two straight elements, a constant-contour, full-span wake using curved 
elements will be more efficient than the straight-line lattice approach since the reduction in 
the number of elements outweighs the increased computation time required by BCVEs. 

As noted above, by choosing each contour filament to be of equal strength the 
filaments are closely spaced in regions of high circulation, and are widely spaced in regions 
where the circulation is low. This relationship between filament spacing and sheet 
strength, and the fact that the filaments follow the actual vortex lines, means that the 
structure and relative strength of the wake can be seen visually. Since the computational 
filaments are aligned along the true vortex lines in the wake, this seems to be the most 
natural and accurate representation of the wake. Also, because the "shed" cross-filaments 
present in the lattice model are eliminated, an additional discretization along the wake length 
is not needed and the potential problems of handling these cross-filaments in a free wake 
representation are avoided. 

Once the wake is laid out, the analysis proceeds much as the tip vortex analysis 
described in Reference 10, except that multiple filaments are trailed from each blade. This 
inclusion of the wake trailing from the full span of the rotor blade gives this method its 
most concise and descriptive names, and moreover is one of the most important features of 
the new model developed here. Subsequent calculations will show the important role 
played by the wake of the inboard portion of the blade. In addition, this full-span approach 
naturally allows for the appearance, evolution, and disappearance of the wake of opposite- 
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sign circulation zones; such zones are very often present near the advancing tip of 
conventional rotors (as well as on highly twisted blades - such as tiltrotor blades - in 
transition flight) and are often critical to the correct handling of wake-induced velocities. 


While offering superior efficiency and resolution of the wake structure, this new 
full-span w ak e approach does introduce certain differences and difficulties. Two new 
features of the wake representation that must be addressed are the formation of loops in the 
wake and the spanwise variation in filament release points as the blade rotates. Handling 
both of these features requires careful treatment of the numerical logic involved in filament 
generation and amalgamation. 

In its current form, the analysis starts with an estimate of the spanwise bound 
circulation distribution at a set of azimuthal locations. This distribution can be calculated 
from a loads analysis, using a simple inflow model, such as the skewed, prescribed 
down wash of the Drees model described in Reference 21 (see Section 6 for a more detailed 
discussion of the initialization and operation of the analysis). Typically, the bound 
circulation is calculated at fixed spanwise stations. The spanwise release points for die 
individual discrete filaments are found by linear interpolation between these fixed spanwise 
calculation points. This must be done at each azimuthal station since the circulation 
distribution changes. Not only do the release points shift along the blade span but some of 
them also disappear and reappear as the maximum circulation decreases and increases. 
This disappearance and reappearance corresponds to the formation of loops in the wake. 

In the initial versions of this forward flight analysis (i.e., prior to the work of Ref. 
14) the bound circulation distribution was assumed to be fixed and did not evolve along 
with the wake geometry. As of the the publication of Reference 15, the restriction on blade 
circulation had been removed, though at that time a simple blade dynamics model was in 
place that precluded fully realistic load calculations. In the current procedure, the free wake 
evolution provides an updated inflow distribution which is used as input to the vortex 
lattice model to generate a new bound circulation distribution. This allows a new set of 
release points to be calculated for further free wake calculations that update the blade 
loading. An interior iteration is required to find a blade dynamics and trim solution that is 
consistent with this blade loading. This will be discussed later, along with details of the 
finite element structural deformation analysis that has been implemented to allow a more 
complete description of the blade motion. 


3.2.1 Wake Structure 

The wake is divided into four distinct zones, as sketched in Figure 3-7 (note that the 
vertical dime nsion is expanded by a factor of 3: 1 in this figure). In Zone 1, filaments trail 
from the inboard ('positive') circulation distribution; no filaments are trailed from any 
inboard reverse flow region. Zone 2 filaments trail from the outboard, negative circulation 
distribution, if one exists at the azimuth in question. Zone 3 is a varying strength vortex 
pair that represents the circulation downstream of Zone 1, beginning at the end of Zone 1 
filame nts and extending a specified azimuthal distance downstream. Zone 4 is the same as 
Zone 3 except that it captures Zone 2 circulation; this zone is rarely important, but is 
included for completeness. 

The filaments in Zone 1 are constant in strength along their length and, moreover, 
the strength of each filament is equal. Varying circulation strength in the Zone 1 wake is 
accounted for by increasing and decreasing the number of filaments present, not by 
increasing and decreasing their strength. For example, if Zone 1 filaments have a 
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Figure 3-7 . Three-view of wake structure trailing from one blade of a 
four-bladed rotor at advance ratio 0.35 . 
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circulation strength of 20 ft 2 /sec and there is a maximum of 100 ft 2 /sec of circulation in the 
bound circulation within Zone 1 at a particular azimuth, then there will be ten Zone 1 

filaments present at this azimuth, five with a strength of -20 ft 2 /sec and five with +20 
ft 2 /sec. Even if the user selected 12 filaments for Zone 1 as input, there will still only be 
10 at this azimuth. The other two will appear at azimuths that have greater circulation. 

At present, a fixed assumption of the wake layout logic is that the radial distribution 
of bound circulation consists of a positively loaded region inboard, with the possibility of a 
negatively lo ade d zone near the tip at some azimuthal stations. Zone 2 filaments will be 
created and deleted in the same manner as Zone 1 for those azimuths at which negative tip 
loading is present. Using this approach for both zones, the filaments represent contours of 
constant vorticity strength in the wake. This procedure could be extended to accommodate 
more than two zones, but it was judged that the present arrangement covered a sufficiently 
broad range of realistic rotor operating conditions for current purposes. 

Unlike Zone 1 and Zone 2 filaments, the strength of Zone 3 and Zone 4 filaments 
vary around the azimuth. Zone 3 contains one varying-strength vortex pair. The strength 
of the Zone 3 filaments at each location is equal to the maximum circulation at that location. 

For the ex amp le above, the two Zone 3 filaments would have strengths of ±100 ft 2 /sec at 
an azimuthal location with a peak circulation of 100 ft 2 /sec. Zone 4 operates in the same 
manner for the negative tip distribution if it is present 

The filament strengths themselves are determined in the following manner. After an 
initial estimate of the blade motion has been made (see below), the circulation on the blade 
as a function of radial position is known at each azimuth. The strength of Zone 1 filaments 
is then set so that they can account for the maximum circulation that will occur at any 

azimuth. For example, if the maximum circulation at any azimuth is 100 ft 2 /sec and the 
user has selected 12 Zone 1 filaments (six will be positive and six negative), then the 
strength of each filament will be ±(1.2)(100)/(6) = ±20 ft 2 /sec. The factor of 1.2 is 
introduced to allow for the possibility that the peak circulation will grow as the calculation 
progresses; lower peak circulations will be automatically accommodated by vanishing 
filaments and by the formation of vortex loops. This 20% "safety margin" has proved 
adequate for nearly all calculations to date. Thus, as the solution evolves, this will be 

enough circulation to account for a peak circulation of up to 120 ft 2 /sec. However, if after 
subsequent evolution of the wake-induced velocity and the blade motion the circulation 
increases above 120 ft 2 /sec then the filament strengths are increased to compensate. 
Convergence of the wake solution is usually retarded (but not precluded) if such a change 
in filament strengths occurs. The same procedure is used to determine the filament 
strengths of Zone 2 filaments based on the (algebraic) minimum negative circulation peak 
outboard of the inboard circulation peak. The 'outboard' stipulation effectively eliminates 
the reverse flow region from the analysis. 

In general, then, the strength of the filaments trailing from Zone 1 will be 


1.2 max(r(r,\|r)) 

1 “ 0.5 Ni f (3.9) 

where Nif is the maximum number of filaments to be trailed from Zone 1, and r(r,y) is the 
bound circulation distribution within this zone. (Note that the factor of 0.5 is added since 
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the total number of filaments that may be trailed from Zone 1 will include an equal number 
of vortices of opposite sign.) The strength of Zone 2 filaments is determined in exactly the 
same way. Unlike the filaments in Zones 1 and 2, the Zone 3 and 4 filaments vary in 
strength along their length. The strength of the Zone 3 filaments at an azimuthal location 

Vo is 


r 3 (vo)= ± max(r(r,Vo)) (3. io) 

That is, the Zone 3 filament circulation is set at the maximum value for the bound 
circulation at the azimuth at which the element was trailed. The circulation strength of the 
Zone 3 trailers will thus vary around the azimuth; the curved elements within this zone are 
assumed to vary in strength linearly over their length. The Zone 4 filament strength is 
defined in an analogous manner. TTie determination of the vortex strength is, of course, 
only one part of the wake layout procedure; the specification of vortex release points is 
addressed in the next section. 

3.2.2 Filament Release Points 

Filament release points are the points on the blade from which the vortex filaments 
are released into the wake to account for the vorticity trailed from the blade. Currently 
filaments are released from the trailing edge of the blade at radial locations that are based on 
the circulation on the blade and the strength of the filaments. (As will be described later, 
special provision is made for a so-called 'overlap near wake' that smooths the flow field 
generated by the vortices in the near wake). Release points vary from azimuth to azimuth 
because of the temporal variation in bound circulation. 

The general formulation for the positioning of release points illustrated can be 
expressed as follows. Denote the circulation distribution along a blade at a particular 

azimuth angle V as r(r), the vortex sheet trailed into the wake has strength -dT(r)/dr. 

Assume for the purposes of discussion that T(r) is positive across the domain of interest, 
so that it corresponds to Zone 1 of the rotor blade; the discussion to follow applies equally 
well to Zone 2, except that the sign of the bound circulation and the circulation in the wake 
are reversed. Assume also that we examine an interval r a < r < rb such that 

r(rb) - no = Ti (3.11) 

where Ti is a specified value of circulation. Tj will be the net circulation around the wake 

trailed from this region. If the interval bounded by r, and rb is sufficiently small, then T(r) 
can be assumed to vary linearly within it. In this case, the centroid of the trailing vorticity 
will be 


f 1 

Rrb) - IXrjJ,, * 2 

Moreover, assuming a linear distribution of T over this interval implies that 


(3.12) 
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r(f) _ rfra) + r(r b ) 


(3.13) 


Thus, if the vortex sheet trailing from the interval r a < r < rb is replaced by a single vortex 
filament of strength Fj, it should be released from the radial location r to ensure that both 
circulation and the first moment of vorticity are conserved; the release point should coincide 
with the radial location of the mean circulation over the interval. 

Figure 3-8 shows a schematic of the release point selection for a representative case 
that illustrates this procedure at work. This case features eight filaments of strength T i 
tr ailin g from Zone 1, while four of strength T 2 trail from Zone 2. (Note that T 2 is not in 
general equal to Ti). The maximum bound circulation here is 3.5Fi. Release points are 
located at the radial locations along the blade where the circulation strength is 0.5Fi, 1.5rj, 

2.5T1, and 3.5Fi. There are eight release points in all, four on the side of increasing 
circulation for the negative strength filament trailers and four on the side of decreasing 

circulation (at positions corresponding to 3.5Ti, 2.5Fi, 1.5Fi, and 0.5Fi) for the positive 
strength filament trailers. The net circulation trailed into the wake is always zero. On the 
scale shown in Figure 3-8, T 2 is half the value of T j, so only two trailers will be released 
on either side of the negative peak circulation in Zone 2. 

In practice, the release point is computed by finding the radial locations at which the 
bound circulation reaches levels that bound the half-integral multiples of Ti (i.e., 0.5T1, 

1 . 5 T 1 , etc.) and then linear interpolation is used to defined the point inside the interval 
where the release is to take place. Consider for example the release of the nth filament from 
the inboard side of the Zone 1 circulation distribution. The bound circulation is in general 
known at a discrete number of radial positions; assume that two stations r a and have been 

identified such that Tfra) < (n + 0.5)r 1 < Tfrb). Then, the nth release point (r v )n is 
dv)n - r, + [ (n+0.5)r 1 - r(r,)] , * Ib ' r *> . 

(Tfo) - Tfra)) (3.14) 

For those special cases at the edges of the Zone, the lower bound is taken to be the point at 
which the bound circulation passes through zero. 

As is noted above, this procedure was developed assuming that the bound 
circulation distribution was piecewise linear, and for good accuracy it is advisable to use a 
sufficient number of filaments along the span so that this condition is approximated 
reasonably well. The consequences of using too few filaments are clear from the simple 

case illustrated in Figure 3-8, since, for example, four trailers with a total of 4.0r 1 of 
circulation strength are trailed from Zone 1, which actually contains a peak circulation of 

3.5^. This discrepancy is clearly related to the use of a relatively small number of trailers, 
and will vanish as the number of trailers is increased. 

Although only a single azimuth station (denoted time step n) is shown in Figure 3- 
8, the figure nonetheless suggests the mechanism by which filament loops are formed. The 
loop in the wake of Zone 1 indicates that at time step n-1 ten vortex filaments were trailing 
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Figure 3-8. Radial bound circulation distribution for a rotor with positive 
and negative circulation zones, showing the release points of 
trailing wake filaments. 

= filament strength in zone 1 
Ty = filament strength in zone 2 












from this portion of the blade, and thus the peak bound circulation was in excess of 4.5Fi. 
The decrease of the peak value to level shown in the figure at time step n has caused two 
release points to disappear and has thus prompted the formation of a filament loop in the 
wake. Conversely, the structure of the wake of Zone 2 in this figure shows that two 
release points appeared at time step n-1, indicating that the negative peak circulation passed 

through -1.5r 2 during this interval, leading to the opening rather than the closing of a loop. 

It is possible that a release point location will appear more than twice along the span 
if the circulation distribution does not increase and then decrease monotonically (i.e., if 
there are local maxima or minima within the zone). The current analysis was tailored for 
the case of a bound circulation distribution with a single positive peak and a possible 
negative peak outboard, and so ignores these variations. Release points are determined by 
moving outboard from the root of Zone 1 for the side of increasing circulation and moving 
inboard from the tip of Zone 1 for the side of decreasing circulation. In this way, no 
release points will be assigned to small 'dips' in the middle of the distribution if they exist. 
The current programming logic could be extended to take account of such features, but this 
was not judged to be necessary in the present effort (though this may ultimately be 
necessary to model fully general load distributions on rotors). Zone 2 release points are 
done in a manner exactly analogous to Zone 1. 

Zone 3 release points do not occur at the blade, but rather at the downstream end of 
Zone 1 filaments. As discussed above, Zone 3 has two filaments of opposite sign whose 
magnitude at any azimuth is set by the maximum circulation strength at that azimuth. All 
the circulation in the wake is rolled up into these two trailing filaments, which are assumed 
to be filaments with linearly varying vortex strength. In rolling up the full-span free wake 
in Zones 1 and 2 to form the varying-strength trailers in Zones 3 and 4, the leading points 
in the rolled-up filaments (sometimes denoted "free wake extensions" in the discussion that 
follows) are placed at the mean of the distribution of vorticity formed by the downstream 
ends of the full-span wake filaments. Since all of the endpoints of the Zone 1 trailers 
should be nearly coplanar, the point representing the mean will lie within this same plane. 
The root filament of Zone 1 is formed from the amalgamation of the filaments trailing from 
the inboard half of this Zone. This Zone 3 root filament will start at a point defined by 
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(3.15) 

Here, x v , y v , and Zy denote the position of the endpoint of the kth trailing vortex from Zone 

1; Yo is the azimuth angle at which the Zone 1 wake terminates. This procedure is repeated 
for filaments indexed (Nif/2 +l)to Nif to form the rolled-up tip filament in Zone 3. Exactly 
analogous procedures are used to define the starting point of the filaments in Zone 4. 

Though this transition to a rolled-up wake can be abrupt in some cases, the flow 
field experienced by the rotor blades will not be adversely affected if the region of the wake 
using the full-span trailing filaments is long enough to move the transition point one to two 
rotor radii from the disk. This is easily accomplished in high speed flight, since the wake 
is in general convected rapidly downstream of the rotor. Low speed cases will feature 
closer interaction of Zone 3 with the rotor, but here the wake filaments from Zone 1 - the 
tip filaments in particular - should already be tightly concentrated by the point where the 
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transition to Zone 3 takes place. The role and formulation of the wake model in Zones 3 
and 4 is discussed further in Section 3.3. 


3.2.3 Wake-Induced Velocity Calculations 

As the blade rotates, the filament points move from their release points due to 
contributions from the convection of the free stream, the velocity induced by die vorticity in 
the wake, and induced velocity from the circulation on the blade. At present, two options 
are available for treating wake-on-wake interactions in the analysis: allowing such 
interactions to include the velocity induced by wakes trailed from other blades or to include 
only the velocity induced by the wake of a given blade on itself. The former is more 
consistent with the spirit of true free wake analysis, and is the preferred approach at low 
and moderate forward speeds where wake distortion is very pronounced. The latter, more 
simplified case is computationally more efficient and is often an adequate approximation for 
high speed cases. 

In either case, the rotor blades experience the velocity field induced by the full rotor 
wake. As will be described in the next section, a vortex lattice analysis is used to compute 
the loads on the rotor blade. The velocity induced by the vortex wake BCVE's at the 
control points within the lattice of vortex quadrilaterals are stored for use in the blade loads 
calculations (see Section 4). The computation of the wake influence at each of the edges of 
the vortex lattice is a computationally demanding task if high lattice density is used, and so 
certain simplifications are invoked, taking advantage of the results of substantial numerical 
experimentation. If the distance from the BCVE midpoint to the quad control point is 
greater than five times distance from the edge midpoint to the control point, then the wake- 
on-blade velocity calculation at that particular edge is bypassed and the velocity at the edge 
midpoint in question is set equal to the velocity at the control point of the quadrilateral. 

A topic of special interest is the computation of the velocity induced by the near 
wake immediately downstream of the generating blade. The basic model used here involves 
the replacement of the curved vortex elements that are trailed from the blade with an 
'overlap region' that produces a smooth flow field at the blade surface. The wake in the 
overlap region consists of extensions of the bound trailers that leave the blade and are 
trailed downstream . The two principal options that exist within the current analysis allow 
the overlap to consist either of straight trailers that are rigidly fixed parallel to the chord line 
of the blade or of trailers that move in accordance with the local free stream velocity. 
Figure 3-9 illustrates a typical configuration of near-wake trailers at four azimuth locations 
for a rotor at advance ratio 0.4. The change in orientation of the trailers around the azimuth 
reflects the varying orientation of the effective local free stream. In addition to the 
alterations in the near wake, a special treatment of the vortex trailers in the overlap region is 
made in the case of highly yawed flow to correct for tip loads, as will be discussed in 
Section 4. The technical substance of the overlap model and the model problems carried 
out to validate it will also be discussed at that point 

In the current overlap model the azimuthal extent of the overlap region may be 
specified. The overlap is only applied to the wake trailing from the generating blade, i.e. 
only points of evaluation on the generating blade 'see' the straight trailers in the overlap 
region. All points of evaluation in the wake or on other blades experience the flow field 
induced by the full free wake. To an observer on the generating blade itself, the free wake 
of filamentary trailers begins at the downstream end of the straight trailers; this has the 
desired effect of smoothing out the discreteness that would otherwise be present if the 
curved filaments began directly at the trailing edge of the blade. 
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Figure 3-9. Schematic of bound vortex lattice and movable near-wake 
trailing vortices in the overlap region: four azimuth angles for 
a rotor at advance ratio 0.4 . 






3.2.4 Vortex Core Modeling and Self-Induction Cutoff Distance Selection 

As with all Lagrangian wake models incorporating vortex filaments, a vortex core is 
required to remove singularities from the flow field. Previous investigations (e.g.. Refs. 
2, 3, and 8) found that load calculations are sensitive to the particular core size chosen. 
One of the objectives in formulating the current model was to remove as much of this 
sensitivity as possible. As long as core radius parameters exist within wake calculations, it 
will be possible to 'dial' or adjust the predicted loads. The discussion that follows is 
directed at outlining the approach that was taken here to implement a treatment of the vortex 
core that constitutes a reasonable step towards reducing the modeling role of vortex 
filament cores. 

It is important to appreciate that the full-span free wake model itself contributes 
substantially to the aim of weakening the modeling role of the core. Alternative models, 
such as using single free tip filaments, must of necessity use adjustments in the core size to 
capture flow field effects which are in fact attributable to spanwise and azim uthal loading 
changes. Since such changes are automatically captured with the current full-span wake 
analysis, one possible ambiguity has been removed. However, since filamentary vortices 
are in fact still used, some effective core structure must be imposed to remove the flow field 
singularities associated with line vortices. The core model used here is the one originally 
proposed by Scully (Ref. 2); the swirl velocity profile is 


v e 


= X. r 
2rc r 2 + r 2 


(3.16) 


This constitutes an 'algebraic' core model, which retains half of the vorticity associated 
with the vortex inside the core radius r c and leaves half outside. Scully notes that this swirl 
velocity profile has considerable experimental substantiation. Also, this particular core 
structure is a convenient choice for the implementation of vortex elements based on ANM, 
as discussed below. 


The issue of the selection of the core radius itself remains. Currently the core r adii 
vary from filament to filament and along filaments from azimuth to azimuth. In keeping 
with the spirit of the discretization of the wake, which places curved filaments on contour 
lines of constant strength of the wake sheet trailing from each blade, the local core radius is 
based on the distance between vortex release points at each azimuth. The core radius 
becomes the average distance to the neighboring release points on either side (at a given 
azimuth) with special cases at the root, tip and center of each zone. That is, for filaments 
numbered n-1, n, and n+1 trailed within Zone 1 or Zone 2, the core size of filament n is 
computed as 


\ — (frv)n * (ty)n-l) ((fy)n+l ' frv)n) 

(chl 2 (3.17) 

At the center, averages are taken between the nearest filament and the maximum circulation 
line. At the root and the tip, averages are taken using the distance to the blade cutout, 
blade tip, or the boundary between Zones 1 and 2. 

In the current analysis, bounds can be placed on the core size and, if desired, 
particular core radii can be chosen for each filament. Numerical experimentation has 
shown that the rotor loading still shows substantial sensitivity if core sizes are adjusted 
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arbitrarily. However, the primary mode of operation for the correlation runs discussed 
below was to allow the core radius assigned to each filament to adjust itself to local 
conditions as outlined above. It is judged that this approach is consistent with the overall 
aim of removing as much arbitrariness as possible from the analysis of rotor airloads. 

In addition to selecting the core properties of the wake filaments, the analysis must 
also compute an appropriate choice for the cutoff distance associated with the self-induction 
velocity generated by each filament. The vortex filament dynamics are currently analyzed 
by integrating the Biot-Savart law over the vortex filaments. This integration is done 
numerically in the present analysis by breaking the filaments up into simpler vortex 
elements (BCVE's in this case) whose integrations can be done in closed form, and 
summing all the element contributions. If a finite-strength vortex filament of infinitesimal 
cross-section is used, then the Biot-Savart integration is logarithmically singular when the 
point of evaluation is placed on the curved filament itself. This logarithmic singularity of 
self-induced velocity is associated with the local curvature of the filament at the point of 
evaluation. This problem is avoided by stopping the integral at a cut-off distance on either 
side of the point of evaluation. This approach is handled by a special self-induction vortex 
element called the SIVE. 

Typically, the choice of cut-off distance for a vortex trailing from a wing or rotor 
blade is obtained from an asymptotic analysis that give results in term of local core size 
(length scale) and the net swirl and kinetic energy content of the vortex core that must 
actually exist at the center of the filament. The energy content and core size can be obtained 
from a detailed knowledge of core structure or, in die case of a tip vortex, by a knowledge 
of the blade load distribution using a recently developed method (Ref. 22). The latter 
approach was used in recent work on the prediction of hovering rotor wake geometry and 
performance (Ref. 23). 

The self-induced velocity for the filaments in the current analysis is also handled by 
a cut-off distance approach, but in this case the cut-off distance is based on the fact that the 
filaments represent a vortex sheet. It can be shown that a curved segment of the vortex 
sheet experiences a self-induced velocity, similar in principle to that experienced by a 
vortex core. A planar, constant strength curved vortex sheet can be shown to have a cut- 
off distance equal to half the sheet width. The velocity computed using the SIVE element 
and this cut-off distance applies at the center of the sheet If the sheet cross-section is 
rotated 90*, as if the sheet now lies on a cylindrical surface, then there is a change in the 
cut-off distance. However, the logarithmic dependence in the SIVE formula makes this 
change insignificant Therefore, in the analysis, the half- width cut-off distance formula is 
used as the default choice for all filaments, though other choices can be imposed if desired. 


3.3 Far Wake Modeling 

The full-span wake described above could in principle be used to discretize the 
entire rotor wake, but this would be computationally inefficient Typically three or more 
f ull revolutions of free wake are required to obtain an adequate model of wake distortion, 
while a full semi-infinite wake is rigorously necessary to capture the complete wake- 
induced velocity field. The high-resolution full-span wake model is required only for the 
regions where wake-on-blade interactions or careful treatment of wake roll-up is necessary. 
(For present purposes, the flow field induced on the rotor itself is the primary concern, 
while in other applications - e.g., the computation of wake-induced flow fields at the tail 
rotor or empennage - a larger region of full-span wake may be required). Beyond such 
regions, a simplified model is appropriate, and the discussion above has given some 
indication of its form. The approach taken is to amalgamate like-sign vorticity trailing from 
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each of the wake filament zones into single, freely distorting tip and root filaments. As 
long as this amalgamation (see Figure 3-7) takes place reasonably far downstream of the 
rotor itself, its effect on the rotor loading will be small. 

Presently, the maximum number of filaments that may be trailed from a given 
azimuth in each of the four wake zones is specified before starting the calculation. This 
number is ordinarily at least ten, but the useful range can be from six to twenty, depending 
on the advance ratio and rotor loading. Smaller numbers of filaments can often be used at 
low speeds, while the complex vorticity field characteristic of high speed flight may require 
a large number of filaments. The analysis also requires the specification of the number of 
turns of full-span wake that are to be used and the azimuthal location at which the two- 
filament wake is to commence. The strength of the rolled-up filaments is adjusted at each 
azimuthal station to conserve circulation in the wake. This accounts for the effect of the 
'shed' component of the vortex wake in the far field. The root and tip filaments in this 
zone (Zone 3, as discussed in Section 3.2.2) originate at radial positions corresponding to 
the centroids of the truncated filaments that are trailed inboard and outboard, respectively, 
of the peak circulation on the rotor blade. These filaments, which will be referred to as 'free 
wake extensions', distort freely as the calculation progresses. 

To complete the analysis of the wake, the free wake extensions must themselves be 
effectively extended in a semi-infinite fashion; simple truncation of the wake can lead to 
errors in induced velocity in some flight conditions, particularly at low speed. The far 
wake model consists of two regions, as shown schematically in Figure 3-10 and is similar 
to the far wake treatment originated in Reference 10. (It should be noted that a precondition 
for invoking this far wake model is the presence of at least two full turns of the two- 
filament free wake extension). The first portion of the far wake consists of a specified 
number of turns of prescribed wake, turns which duplicate the geometry of the last turn of 
the free wake extensions. The spacing and orientation of these prescribed turns is set by 
the average convection over the last two free wake extension turns. The prescribed turns 
consist of BCVEs, and since their geometry depends on the configuration of the last turns 
of free wake, this prescribed wake evolves constantly during the free wake evolution. 

Given this model for the far wake, the problems to be resolved are the 
determination of the direction in which the wake is convected, and the development of an 
efficient method to sum the infinite number of repeated turns of far wake. The mean 
convection direction is determined by comparing the positions of the last and next to the last 
turns of free wake. The change in position of each set of corresponding points on these 
two turns is computed. These position differences between pairs of points, originally 
generated at the same blade azimuthal location, defines a displacement vector. A running 
average of all these displacement vectors, generated while stepping the blade through the 
previous full turn, is computed to determine an average displacement vector. This average 

displacement vector, r s , is taken to represent the average displacement per blade 
revolution of a wake turn due to overall convection. It is used to orient the successive far 
wake turns in space. The purpose of the running average over a cycle is to filter out the 
distortion effects that are superimposed on this mean convection. 

Note again the far wake configuration illustrated in Figure 3-10. The far wake can 
be thought of as being composed of a set of semi-infinite rows of elements running in the 

direction of U , one of which is shown in the figure. It is possible to sum the 
contributions of each row in an efficient manner by using a discrete approximation for the 
Biot-Savart law. )This development is discussed in detail in Reference 10 and will not be 
repeated here.) 
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Point of Evaluation 



Figure 3-10. Schematic of the far wake model and its coupling to the 
two- filament free wake extensions. 
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This intermediate prescribed wake model is based on die idea that the best estimate 
of the shape of a far wake turn is the shape of the turn of the free wake extensions. Ideally, 
after each blade revolution, the vortex wake represented by the last turn of the extensions 
would occupy the position of the first far wake turn. This assumption corresponds to 
freezing the shape after the last turn of the extensions. The effect of subsequent 
downstream distortion is, therefore, not included, and only convection in an average sense 
is taken into account. This approximate shape is best for the first few far wake turns, since 
there has been relatively little time for subsequent distortion. Neglecting this distortion in 
the far wake turns that are farther downstream makes little difference since the influence of 
their detailed shape is less crucial because they are farther away. 

The far wake model has been implemented into the RotorCRAFT analysis and has 
performed in a fashion similar to previous forward flight wake studies (e.g., Ref. 10 and 
1 1). A smooth, stable transition between the end of the free wake and the beginning of the 
far wake has been achieved. Since the far wake structure is based on free wake 
information, the far wake shape evolves along with the free wake solution. Note that there 
is a different far wake shape at each time step in the solution. The far wake configuration 
also adjusts automatically to the rotor flight condition. 

3.4 Implementation of Wake Models Based on ANM 

In the free wake analysis of rotors in forward flight or hover, the use of BCVEs has 
led to considerable success, as is well documented in the literature ( Refs. 9-15). The 
proven robustness and accuracy of the method are direct consequences of its high 
resolution representation of the filamentary wake vorticity field. However, in some 
applications, the use of such high resolution, computationally intensive vortex elements 
everywhere in the flow domain may significantly reduce the efficiency of the method; while 
a high resolution representation of the vorticity field is required to capture the steep velocity 
gradient when an evaluation point is in the immediate vicinity of a vortex element, the bulk 
of the calculations usually involve evaluation points which are in the far field where coarse 
resolution is sufficient. In this regard, the use of the particle-based "vorton" method can 
result in considerable savings in computational costs. Because vortons are "clumps" of 
vorticity concentrated at a single point in space, the algebraic functions associated with the 
calculation of their velocity field are particularly simple and inexpensive to evaluate. 

In the literature, applications of this method to unsteady vortical flow problems 
have shown promising results (Refs. 24-26). However, convergence studies of the method 
(Refs. 27-29) indicate that sufficient overlap of the particles must be maintained, i.e. 

5s = O , where & is the spacing and O is the core size. This requires that a large 
number of vorton particles be used. For a typical rotorcraft free wake calculation where 
convective washout of an assumed wake shape is used to obtain converged solution, the 

order of the calculation is O(N^) , where N is the number of vortex elements, and for a 
large number of vorton particles, this represents a formidable computational load even with 
the simple functional form of the individual velocity calculation. 

Recently, a new method which incorporates both the accuracy of BCVEs and the 
efficiency of a vorton method has been developed (Refs. 16 and 17). In this method, 
referred to as the method of Analytical/Numerical Matching (ANM) , a low resolution far 
field solution is constructed with a small number of fat vorton particles. The small number 
of elements limits the computational load to within the capabilities of conventional 
computing machines while the fat core size assures the mathematical consistency of the 
scheme. For cases of close approach, where an observation point is in the near field of a 
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vortex filament, a local analytical correction is applied to capture the high gradient near field 
solution (Fig. 3-11). Inspired by the method of matched asymptotic expansions, the local 
analytical correction serves to match the coarse resolution far field solution to the fine 
resolution near field solution and is composed of two opposite-signed curved arc elements 
which are fitted locally to the vortex filament. The positive sign arc has the actual physical 
core size while the negative arc has the fat core and is used to match the fat core particle 
solution with the thin core arc solution. In principle, any curved arc element can be used for 
the local correction; for present purposes, however, a closed circular ring is chosen because 
of its closed-form Biot-Savart integral which ensures the efficiency of the velocity 
evaluation. Furthermore, the correction is applied only as needed for near field interactions, 
which usually represent a small fraction of the calculation. The resultant scheme is highly 
efficient, providing the accuracy of the BCVE method but with the CPU savings of a 
vorton method. 


3.4.1 Results of Test Runs 

In their work. Bliss and Miller (Ref. 17) have applied the ANM method to several 
test problems involving simple geometry such as circular and elliptical vortex rings. They 
were able to obtain good agreement between the calculations and theory with an impressive 
speed-up factor of three to four compared to corresponding BCVE calculations. In an effort 
to independently establish the robustness and accuracy, as well as the parameter sensitivity 
of the ANM method, extensive test runs using the ANM method have been carried out as 
part of the development of the RotorCRAFT analysis. The parameters associated with the 

method are 5s , Gf and &n. , where Ss is the spacing between vorton particles, is 
the fat core radius and OJi is the near field cutoff distance, within which the analytical 
correction is applied. The following is a brief description of some of the different test 
cases. 


The first sample calculation involved a circular vortex ring with circulation i"= 1 , 
radius R= 1 and core size 0=0.2 discretized into 50 equally spaced vortex elements. The 
induced velocity normal to the plane of the vortex ring was computed at evaluation points 
along a radial line in the plane of the ring. Both ANM and BCVE (with the self-induced 
SFVE component included) formulations were used and the results have been compared. 
The ANM parameters used are Gn=5Gf , &f=Ss and &=2 teR/ 50 . Figure 3- 12a shows 
a comparison of the normal velocity computed along a radial line that intersects a vortex 
particle. Very good agreement between the two calculations can be observed. A plot of the 
percentage velocity error, given in Figure 3-12b, shows that the error of ANM relative to 

BCVE is essentially zero in the region 0.4 <R< 1.6 , where the analytical correction is 
applied. A slight discrepancy is observed at /?=1.0 , the center of a vortex particle, and 

this is due to the difference in the self-induction model. In the region R^OA or R> 1.6, 
slight discrepancies which are within an acceptable level, are observed. In these 
calculations, a speed-up factor of four for the ANM calculation was realized. 

Additional calculations were run for a 4:1 elliptical vortex ring with 7=1 and 
0=0.2 discretized into 50 elements. Cosine spacing is used to give a high density of 
elements in the region with steep curvature. The velocity normal to the plane of the ring is 
computed on the major axis, where the error of the ANM method, if any, is expected to be 

a maximum. Figure 3- 13a shows a comparison between the ANM calculation, with Gf=8s 
and Qn=50) , and the BCVE calculation. Significant over-prediction of the velocity near 
the interior of the ellipse, as reported in Reference 17, is observed. Calculation of the local 
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igure 3- 11 a. Three-part ANM composite solution: (1) fat core vortex 
particle outer solution, (2) thin core analytical inner solution, 
and (3) opposite sign fat core analytical matching solution. 
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Figure 3-1 lb. Schematic of ANM composite solution. 
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radius of curvature shows that it is less than the fat core used, implying that the fat core 
analytical matching is given by a small disc of vorticity, i.e. a vortex ring with core size 
greater than its radius. Clearly, this is not a suitable model for matching the fat core particle 
solution to the local thin ring solution. Repeating the ANM calculation with a reduced fat 
core of Cfy=0.5& while keeping the near field cutoff parameter unchanged yields 
im pr oved results, as shown in Figure 3- 13b. This suggests that in using the ANM method, 
care must be exercised to ensure that the fat core size used does not exceed the local radius 
of curvature and this can be achieved either by using sufficient density of vortex elements 
in the region of high curvature, or by limiting the fat core size used. Test runs involving 
variation of On have shown little effect, verifying the insensitivity of the results to On , 
which is an advantage for the method. 

As these sample problems indicate, the application of vortex wake models using 
ANM holds considerable promise, in that flow fields with accuracy comparable to BCVE 
calculations can be achieved with considerably reduced computation time. A variety of 
calculations were undertaken to explore the suitability of ANM elements for routine use 
within the RotorCRAFT analysis. The test cases discussed here involve computations on a 
representative four-bladed rotor at both low and high advance ratio. The particular 
configuration used is an H-34 main rotor, operating at a thrust coefficient of 0.0037. 

In the case of low advance ratio, the wake vortex filaments tend to remain close to 
the rotor for a substantial period of time because of the low convective flow speed. This 
results in significant close interactions between filaments and represents a severe test for the 
robustness of the ANM-based wake model in terms of its capability to handle close 
encounters. In this case, one revolution of the H-34 rotor at an advance ratio of 0.15 is 

computed with 16 steps per revolution. The ANM parameters used are 0f=0.5& , 
<7 n =3 <7/ . Figure 3- 14a shows a comparison of the time history of the computed airload at 
the 90% radial station using the ANM and BCVE-based calculations. Comparing these 
cases, it is observed that the ANM calculation tracks its BCVE counterpart closely. 
However, a 50% saving in CPU is realized when comparing the ANM to the BCVE 
calculation. 

In the next case, calculations are carried out at an advance ratio of 0.39 and 
comparison of the airload time history shows excellent agreement (Fig. 3- 14b). For the 
higher harmonics, however, the comparison, as shown in Figure 3- 15a, shows some 
differences in amplitude and phasing. In terms of amplitude, there are significant 
underpredictions by the ANM calculation at the first, second and last minima, and the first 
and last ma xima. There are also phase errors at the second maxima and minima. This can 
be attributed in part to errors that can be introduced at high forward speed due to the 
increase in the length of some of the wake elements; since the size of the fat core on the 
vorton elements is set at one half of the element spacing, stretching of the elements may 
lead to unanticipated increases in the fat core size. However, repeating the ANM 

calculation using a reduced fat core of O/=0.25& improves the amplitude comparison 
significantly, as shown in Figure 3- 15b. Some phase difference remains, but this may be 
due to the coarse time discretization, i.e. number of steps per revolution. Running the 
calculation at a finer resolution of 24 steps per revolution improves this slightly. 


3.4.2 Discussion of Parameter Selection for ANM 

In these sample airload calculations of a rotor in forward flight, convective washout 
of an initially assumed wake geometry is used to step the solution towards convergence. 
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Figure 3-14. Comparison of rime history of airloads at r/R = 0.9 for an 
H-34 rotor using ANM and BCVE wake models. 
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a) ANM fat core equals one-half particle spacing. 



b) ANM fat core equals one-quarter particle spacing. 
Figure 3-15. Comparison of time history of higher harmonic airloads at 
r/R = 0.9 for an H-34 rotor at advance ratio 0.39 using 
ANM and BCVE wake models. 
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For N vortex elements, this usually involves N steps to converge, resulting in an 0(N ) 
calculation. Furthermore, for the blade dynamics calculation to converge, several rotor 
revolutions are required and this represents a significant computational load. To reduce the 
computational time, the ANM method seems a favorable candidate, in that its accuracy is 
comparable to BCVEs while it has a distinct advantage in efficiency. In this effort, the 
ANM method has been added to the forward flight code as an option to reduce 
computational load. Representative test runs, presented in the last section, have shown that 
a speed-up factor of two can be realized. 

Despite the attractive savings in CPU, the ANM scheme should not at present be 
thought of as a replacement for the BCVE scheme, but rather as an alternative in cases 
where a significant time constraint is encountered. This is because as demonstrated in the 
last section, care must be exercised in using the ANM scheme. This is especially true for 
cases of high advance ratio, where wake vortex element spacings are large due to the high 
convective flow speed. In these cases, the fat core size of the particle should only be a 

small fraction of the spacing and test runs have indicate that a fat core not exceeding 0.25 & 
would be appropriate. This represents a sacrifice in the smoothness of the far field solution 
but is essential in assuring the correctness of the near field matching. For cases involving 
low to moderate advance ratio, the significance of such differences are diminished, and a 

fat core of OM).5& should be used to ensure the smoothness of the wake vorticity field. 
In terms of the near field cutoff parameter Gn , test runs have showed little dependency of 
the results on the parameter and in general, a cutoff of G„ =3 Gf will ensure a proper 
matching of the far field and near field solutions. 

In conclusion, it is anticipated that as greater familiarity with ANM wake elements 
is achieved, they will become the computational tool of choice, replacing the existing 
BCVE formulation. For this reason, they have been included as a computational option in 
the current RotorCRAFT analysis, though for the present the selection of BCVEs remains 
the preferred approach. 
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4. VORTEX LATTICE MODELING OF ROTOR BLADE AERODYNAMICS 


As noted above, one of the primary tasks of the Phase I effort described in 
Reference 14 was to introduce a lifting surface/vortex lattice model of the rotor blade into 
the computation of rotor airloads. The principal motivation for this effort was to allow 
improved resolution of loads in the vicinity of the tip as well as during close blade vortex 
interactions (though the current formulation is still inadequate for prediction of loads due to 
direct vortex impingement on the blade surface). This model has been revised and 
extended during Phase II to include a refined treatment of tip effects as well as of both 
steady and time- varying near wake effects. After the initial effort described in Reference 
14, though, a variety of areas for improvements in the baseline model were identified. 
Among these were the need to incorporate profile drag into performance predictions, 
compressibility modeling, the inclusion of lift stall , the influence of yawed flow and wake 
skew on load predictions, and a model of the inflow contributed by the extreme near wake. 
The following sections will discuss the steps that were taken to include all of these features; 
the discussion will overlap at many points since several of the topics are closely related. 
These topics will be preceded by discussion of the background on vortex lattice treatments. 


4.1. Background on Vortex Lattice Models 

Many modem rotors feature complex planforms, and it is advantageous to use 
lifting surface routines (which are a subset of the general class of aerodynamic panel 
methods) to analyze such designs. Subsonic panel methods have been widely used in the 
computation of flow around both lifting and nonlifting configurations (e.g., Refs. 30 and 
31). These methods use a distribution of sources, doublets or vortices, or different 
combinations of these singularities arranged in a variety of geometrical forms. A wide 
variety of panel methods are in routine use in the fixed-wing aerodynamics community, and 
they have been applied to the full range of complex configurations including surfaces with 
thickness, complex configurations involving separated flows, surfaces interacting with 
nonlinear vortical wakes, and the like. 

The particular model employed here focuses primarily on thin lifting surfaces with 
no side- or leading-edge separation (with the exception of tip effects in yawed flow 
discussed in Section 4.7), and so is in the tradition of models first developed by Falkner 
(Ref. 32) and later popularized by such researchers as Rubbert (Ref. 33) and Margason and 
T -amar (Ref. 34). The rotor blade is represented in the present work by a surface consisting 
of vortex quadrilaterals. Four constant strength straight-line vortices form the sides of each 
q uadrilater al, except at the trailing edge of the blade where a modified vortex quadrilateral 
without the downstream line vortex is used. Also, as noted above, special modifications 
can be made to the side edge of the tip in cases where yawed flow causes part of the wake 
to be trailed from this edge. 

Lifting surface theory has received relatively little use in rotorcraft applications, 
with the exception of such work as References 35 and 36. Lifting line theory has been 
preferred for routine use because of its simplicity and because its form simplifies the direct 
use of two-dimensional airfoil data to account for transonic and viscous effects. On the 
other hand, lifting surface theory has the advantage of a more correct treatment of three- 
dimensional effects. Three-dimensionality is important in the localized effects of blade 
vortex interaction and for the prediction of tip loading, particularly for treating novel tip 
shapes with considerable sweep and taper. One objective of the present effort was to 
explore the development of a hybrid scheme in which the advantages of both methods are 
incorporated. 
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4.2. Outline of the Vortex Lattice Layout Procedure 

The lifting surface/vortex lattice analysis used here was adapted for use from the 
hoVer simulation described in Reference 23. This formulation allows substantial flexibility 
in the specification of the blade's planform so that complex designs may be accommodated. 
Currently, the lattice can be divided into as many as ten different regions, with separate 
linear distributions of twist, taper, and sweep within each. The spacing of the 
quadrilaterals in a vortex lattice analysis is an important consideration, as is discussed in 
Reference 37. The judicious selection of the density, spacing, and orientation of the 
quadrilaterals can considerably enhance the efficiency and rate of convergence of the blade 
loading. The current analysis has been provided with sufficient flexibility to arrange 
essentially arbitrary chordwise and spanwise distributions of lattice elements though the 
control points are always assumed to lie at the geometric center of the quadrilateral. For 
simplicity and convenience, however, all of the calculations discussed in this report feature 
uniform spacing of vortex quadrilaterals both in the chordwise and spanwise directions, 
unless otherwise noted. 

The vortex quadrilateral lattice is drawn in blade coordinates, as defined in Figure 
4-1. First the blade segments are laid out separately in the XY-plane applying taper and 
sweep. The lattice is displaced toward the trailing edge by a distance of one quarter of the 
chordwise length of the leading edge quadrilaterals. For one row of quads and an unswept 
rectangular planform, this puts the quad leading edge (load vortex) along the quarter chord 
line of the blade and the vortex lattice control points (center points of each quad) along the 
3/4 chord line of the blade. The lattice is inset from the blade root and tip by a distance 
equal to a quarter of the width of the last quad at either edge. This technique was suggested 
in Reference 38, and was found both there and in sample calculations carried out for this 
effort to accelerate convergence of the loading solution. 

For reference purposes, the quarter-chord line of the blade is taken as the line that 
connects the quarter chord points of each blade section, while the X axis of the blade 
coordinate frame is the line connecting the hub with the quarter-chord of the root section. 
The sweep angle for any segment is defined as the angle the local quarter-chord line makes 
with the X axis. Pitching moment calculations use the local quarter-chord line as a 
reference axis, however both collective and cyclic pitch are applied about the X axis. 

In cases more complex than untwisted rectangular planforms, additional steps must 
be carried out Figure 4-2 depicts the individual steps carried out in laying out the blade 
geometry. In actual calculations, the order of operations is different than is shown here; 
first taper is applied linearly from root to tip along each segment Then sweep is applied 
by displacing each segment toward its trailing edge (+Y-direction); the sweep angle is the 
angle between the X-axis and the quarter chord line. The twist gradient is applied by 
rotating each chord of the lattice about its quarter chord point Finally, anhedral is applied 
by displacing each segment downward in the +Z-direction. The resulting vortex lattice 
structure is stored and written to a file to be used in a graphical verification of the 
planform. 

In addition to the geometric inputs just described, a camber distribution for the 
blade may also be specified. Numerical input is used to describe the geometry of the 
camber line. When camber is present, the lattice itself is not deformed to fit the specified 
distribution, rather the boundary conditions at the vortex quadrilateral control points are 
altered to introduce the surface slope into the calculation. However, a certain minimum 
chordwise density of quadrilaterals is required to resolve the camber distribution; an 
absolute minimum of three quadrilaterals chordwise should be used, with five or more 
being desirable. This level of quadrilateral density can create a substantial computational 
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Blade Coordinates (X,Y,Z) 



Figure 4-1. Definition of shaft coordinates and blade coordinates. 
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burden. An alternative approach for including zero lift angles is also available for cases 
where improved computational efficiency is desired. An appropriate rotation of the vector 
normal to the blade at the control points of each quadrilateral will alter the effective angle of 
attack of the section and can be used to introduce the shift of the zero lift angle of attack. In 
this manner, realistic zero lift angles of attack can be introduced into the calculations with 
only one quadrilateral chordwise. The zero lift angle is computed from two-dimensional 
airfoil data read in to support the computation of profile power, a task described in Section 
4.4 below. That section will also describe the methods used for incorporating pitching 
moment information into calculations where camber distributions on the blade are not read 
in. 


4.3 Computation of Aerodynamic Loads 

The discussion above outlines the procedure used to lay out the vortex lattice to 
represent a general rotor blade planform. The solution method used to find the bound 
circulation given this lattice is essentially a straightforward implementation of the classical 
approach described in the literature on lattice methods for fixed wing applications (e.g., 
Ref. 34). Each of the quadrilaterals is examined individually and a mean vector normal to 
the quadrilateral surface is established as shown in Figure 4-3, which also shows the 
location of the 'control point' associated with the quadrilateral. Given this and the location 
and orientation of each of the quadrilaterals on the blade, the velocity induced by the blade 
lattice on each of the control points is determined, assuming unit strength for each 
quadrilateral. Then the resulting velocity is resolved in the normal direction at each control 
point, yielding an array of influence coefficients relating the vector of bound circulations, 

y , to the downwash, w , at each control point: 


w = Ay 


where 


(4.1) 


A 



i , j = l,....n 


(4.2) 


Here, n is the number of vortex quadrilaterals on the blade. Note that the location of each 
of the quadr ilaterals as well as the velocities above are defined in the blade reference system 
depicted in Figure 4- 1 . 

This array of coefficients is stored and the velocity induced at each control point by 
the free stream, the blade rotation and deflection (discussed below), and the rotor wake are 
summed and resolved into 'normal- wash' velocities that are then used to find the vector y 
of bound circulation values on the disk as follows: 


flj " (<1 free stream) j + (*lwake)j + (<I blade) j 


(4.3) 


w 


j 


= q j n j 


(4.4) 
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Normal vector at X c is n = R j x R 2 


X c = 0.25 (Xi + X 2 + X 3 + X 4 ) 


Figure 4-3. Typical vortex quadrilateral, showing the comer indices and 
diagonal vectors. Control point found using the mean of the 
comers. Normal vector defined as indicated. 
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= A' 1 w 


(4.5) 


The vector y can then be used to solve for the forces and moments on each segment of 
the blade by applying the Joukowsky Law to each of the four edges of the vortex 
quadrilaterals, i.e., 

Fjk = P Yj (qjk x sjk) Aj k , j = 1,—n k = 1 , 2, 3 , 4 (4 6) 


Here, s jk is the unit vector directed along edge k of quadrilateral j ; £jk 

length of this side, while Yj is strength of the quad. The reference velocity q for the 

evaluation of forces is computed at the mid point of the edge k . This velocity contains all 
the components deriving from the free stream, the wake, the motion of the blade, and the 
velocity induced by the entire vortex quadrilateral grid, though the influence of the segment 
on which a particular point lies is deleted. These individual forces are then summed to yield 
the integrated forces on each blade: 

F = (HI S + YJ S + TK S ) = XX F jk 

j k (4.7) 

Moments exerted by the blade about the sectional quarter-chord reference axis can also be 
computed: 

M = (LI S + MJ S + NK S ) = XX F jkxr jk 

j k (4.8) 


Moments about the flapping and lagging axes of the blade are taken to compute the gross 
blade motion, while pitching moments can be computed for each section to use for the 
calculation of torsional deformation. Here, r is the vector from the reference axis to the 
point of action of F. (Tabulated information on pitching moment coefficients can also be 
used to determine the moments, as will be described in the next section). 

The resulting forces and moments are computed in the blade reference frame and are 
then transformed into the shaft frame to yield the thrust and induced torque. Additional 
computations are undertaken to resolve those components that contribute to the horizontal 
and side forces (H and Y, respectively). Of course, these calculations omit the effect of 
viscous forces on the blades, and steps must be taken to include them. However, the basic 
vortex lattice analysis can be modified to include the effect of compressibility on the 
sectional lift and moment, and the approach taken to these issues is described separately 
below. Other modifications are necessary to account for stall, yawed flow, and unsteady 
effects in the near wake; these topics are also taken up in the following sections. 


4.4 Prediction of Sectional Profile Drag and Profile Torque 

To introduce forces generated by profile (viscous and pressure) drag into the 
calculation, the only practical approach at present is the use of two-dimensional airfoil data. 
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The form of data required for this coupling is tabulated values for sectional drag 
coefficients as a function of both Mach number and sectional angle of attack. The definition 
of the local Mach number at a particular radial station is relatively straightforward, being 


M(r,\|r) = C 1 K 


((£) 


|isin \j/cos a s J 




( 4 . 9 ) 


though adjustments for tip relief are required in general. A more uncertain issue is the 
definition of an effective angle of attack for use in the look up tables. 

The definition of effective sectional angle of attack is intrinsic to lifting line theory, 
being defined as 


etc = (X g - cq 


( 4 . 10 ) 


where 04 is the inclination of the effective free stream to the local zero lift line. The 

definition of the 04 in the context of rotary wing lifting line theory is relatively 
straightforward, being 


cq = tan 



( 4 . 11 ) 


where Up is the flow field component normal to the section in blade coordinates, while 
Uy is the tangential flow component. Up itself consists of contributions from the free 
stream, the blade motion, and the wake-induced velocity field. In the context of lifting line 
theory, the wake-induced velocity is evaluated at the quarter-chord point of the section. 
However, in vortex lattice models, there is no clearly defined reference point at which the 
wake-induced velocity can be computed to yield its contribution to the effective angle of 
attack for the section. 

To define such a reference velocity in the context of lattice methods, first recall that 
in the limiting case of a high aspect ratio wing or rotor blade, a good approximation to the 
two-dimensional flat plate solution for a lifting airfoil can be recovered by using a large 
number of chordwise vortex quadrilaterals. In this sense, the portions of the lattice running 
parallel to the blade span constitute the "two dimensional part" of the lattice solution. 
Conversely, the trailing legs of the vortex quadrilaterals represent the "three dimensional 
part" of the solution, since they would all have zero strength on a wing or blade of infinite 
aspect ratio. The flow components associated with these trailers are directly related to finite 
span effects. On this basis it is appropriate to define an effective angle of attack for a 
section by deleting the flow field contributions from the span wise quadrilateral edges and 
computing simply the contributions from the chordwise or trailing legs of the lattice. 

This reference velocity, suitably averaged over the chord for cases involving more 
than one quadrilateral chordwise, provides the component of induced velocity that is 
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contributed by the blade and near wake trailers in the overlap region to the complete 
induced angle of attack. (Note that this approach reduces to the lifting line definition of 
induced velocity for the case of one chordwise quadrilateral). The resulting velocity is then 
added to the velocity field generated by the rest of the wake in the manner indicated in 
Equation 4.10 to yield the effective angle of attack. The local Mach number and angle of 
attack can then be used to enter two-dimensional look-up tables to compute drag and 
moment coefficients; lift coefficients are also routinely looked up, though they are used 
only to define post-stall characteristics, as described below. 

The two-dimensional coefficients are, of course, defined only for a specific airfoil 
section. Many rotor blades feature more than one section along the span. In the current 
analysis, as many as ten sections along the span may be specified. For any given span wise 
station, the section coefficients are computed for each of the two airfoils that bound the 
segment containing the station of interest and then interpolated linearly to the desired point. 
Similarly, for each Mach number/angle of attack pair, bilinear interpolation is used to find 
the appropriate coefficients within each look-up table of lift, drag, and moment coefficients. 


4.5 Compressibility 

Compressibility has an important effect on rotor performance at typical tip Mach 
numbers for modem rotors. Some of this effect is captured by the inclusion of Mach 
number dependence in the look-up tables used for profile drag coefficients. However, 
compressibility also has a significant impact on the lift generated by airfoils at specified 
angles of attack, and so its influence on thrust and induced power must be considered as 
well. This section discusses how the augmentation of lift associated with compressible 
flow can be achieved. Preliminary work has also indicated that some of the nonlinear 
aspects of airfoil section behavior associated with transonic flow (i.e., alterations in 
moment coefficients and zero lift angle) can be incorporated in the context of nominally 
linear vortex lattice methods; implementation of methods to accomplish this has been 
reserved for follow-on work. 

For two-dimensional thin airfoils, treatments like the Prandtl-Glauert correction to 
lift curve slope are acceptable methods for including compressibility effects. In vortex 
lattice calculations, transformations that are similar in spirit but more elaborate in detail 
must be invoked. When analyzing hovering rotors, Kocurek, in Reference 35, used a 
transformation of the entire space surrounding the rotor in generating a correction for 
compressibility. Currently, a more restricted transformation of the blade geometry is used, 
one based on the local Mach number at the radial stations along the span. 

It is assumed that the flow around a given blade section can be found using the 
compressible potential equation 


P 2 *xx + 4yy + *22 - 0 P 2 = 1 - M 2 (r, V ) (4.12) 


where M represents the local Mach number as defined in Eq. 4.8. (For this discussion, 
x denotes the streamwise flow direction (positive downstream); y lies along the radius 
(positive out); and z is positive up.) Using a transformation whereby x is replaced by 

(3xj , y by yj ,and z by Zj changes the above equation to: 
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^xiXj + ^yjyi + ^zjzi " 0 


(4.13) 


This equation governs the incompressible flow around a transformed blade in which the 

blade chord is stretched by a factor of (note that f) varies along the radius). The 
transformed blade is used in the time stepping analysis to find the wake geometry and the 
associated loads on the blade. Then, the loads must be corrected in accordance with the 
geometrical transformation to obtain the thrust and torque on the rotor in compressible 
flow. 


The transformation outlined above leaves the bound circulation, the downwash on 
the blade, and the dimensional lift per unit span the same for both the compressible and 
incompressible problems. Thus, the thrust and induced power computed for the 
transformed, incompressible problems are the same for the original, compressible problem. 
However, because of the stretching of the airfoil chord, the lift coefficient for any section 

for the compressible problem is greater by a factor of than the lift coefficient in the 
transformed, incompressible problem. 


4.6 Stall 

To this point, the discussion has focused implicidy on the prediction of the loading 
at low to moderate angles of attack. In many forward flight conditions, sections of the 
blade span can reach or exceed the static stall angle. A vortex lattice calculation will 
compute unrealistically high sectional lift for such cases, since no allowance is made for lift 
limitation or indeed for other than linear lift characteristics. It is inappropriate to expect that 
models currently available will be able to offer detailed insight into the loading on those 
portions of the rotor blade undergoing stall. However, reasonable steps can be taken to 
ensure that the major features of the lift behavior of rotor blades can be captured. 

First, it should be noted that the presence of yawed flow will substantially increase 
the nominal maximum lift coefficient that a section of the blade may carry. Experimental 
evidence of this phenomenon has been observed in experimental work on fixed wings, and 
the qualitative model used to allow for this effect is discussed immediately below. Also, 
many rotor configurations have a substantial sensitivity to rate effects in delaying the onset 
of stall. Reference 5 among many other sources discusses the potential importance of 
dynamic stall. No dynamic stall model is used in the present version of RotorCRAFT, 
though it is a prominent candidate for addition in follow-on development 

As discussed in Section 4.4, the local effective angle of attack is computed at each 
radial station on the blade, as is the section lift coefficient and the maximum lift coefficient 
permitted by the look-up tables. Once the computed c^ (appropriately adjusted for yawed 

flow) exceeds the stipulated maximum lift coefficient, the stall model is invoked and the 
effective angle of attack is used to look up an appropriate lift coefficient given the Mach 
number of the onset flow. The vortex lattice lift calculation for the section is overridden by 
a new value for sectional lift with the form 


e. =1 p (uf. + U^cc* 
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(4.14) 



This new lift calculation is then resolved to provide thrust and torque contributions in the 
shaft axis system for subsequent use in performance and dynamics calculations. The local 
bound circulation is also reduced to properly reflect the reduced lift at the section. The 
computation of profile drag is not affected, since 2-D look-up tables using the effective 
angle of attack are used for both stalled and unstalled sections. Sample calculations using 
this procedure in the context of wing calculations are discussed in Section 4.9. 

4.7 Overlap Near Wake 

As described in Section 3, the rotor aerodynamics analysis utilizes a free wake 
having constant strength vorticity contours coupled to a vortex quadrilateral model for the 
blade aerodynamics. To achieve consistent results, an overlap wake model is used behind 
each blade to capture the effects of the extreme near wake. The overlap wake is a finite 
length prescribed wake consistent with the vortex quadrilateral model, namely it consists of 
a set of vortex filaments trailing behind the quadrilateral structure. When calculating 
downwash at points on the blade, the overlap wake provides a near field contribution, 
while contributions from the remainder of the wake are provided by the constant strength 
contour free wake. On the other hand, when calculating velocities in the wake, the free 
wake is extended right up to the blade, hence the overlap. The overlap wake is needed 
because the free wake model is not consistent with the blade vortex quadrilateral model in 
terms of either the number of filaments or their location with respect to the blade control 
points. 


The following technical issues were examined with respect to the overlap wake 
model and its role in die calculation of the blade aerodynamic loads: 

Accuracy of the truncation of the overlap wake and its replacement with the constant 
vorticity contour wake beyond a certain distance. 

Treatment of the overlap wake with respect to skewed inflow, i.e. effective sweep, 
as the blade rotates. 

Treatment of the blade tip as downstream edge in the presence of skewed inflow 
(effective sweep). 

These issues were studied by examination of the simpler problem of a straight rectangular 
wing in a skewed (swept) inflow. This problem contains all the physics essential to the 
above issues, without the complexity of a rotor calculation. The important physical effects 
are thereby more easily separated and examined. A simplified version of the same vortex 
quadrilateral code used for the rotor analysis was used for this study of corresponding 
effects on a vortex quadrilateral wing. The study of wake truncation is presented here, 
while the results of the skewed inflow analysis are presented in the next section. 

The effect of wake truncation was examined by considering a rectangular vortex 
quadrilateral wing in uniform, unswept flow. At a certain distance behind the wing the 
tr ailin g vortices were truncated and replaced by a vortex pair. Each vortex in the pair has 
the samp, net circulation as the wing bound circulation, and each is located at the spanwise 
centroid of the corresponding trailed vorticity from the near wake. This case is an extreme 
one, in that it corresponds to a rotor blade with a constant vorticity contour free wake 
having only two filaments. Note that the positioning of the free wake filaments in the rotor 
code is essentially equivalent to centroidal positioning. The geometric layout is shown in 
Figure 4-4 ("centroid model"). 
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Figure 4-4. Schematic of the model problem for the study of wake 
truncation. 



Figure 4-5. Predictions of induced drag on a wing with aspect ratio 4.0 
as a function of y t . 
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The planar vortex wake in Figure 4-4 is composed of N discrete parallel vortex 
filaments, truncated a distance yt behind the trailing edge. Such vortex systems are 
commonly used in the vortex lattice method, or VLM, to model the trading wake of a finite 
wing, except they normally extend to infinity. The vortex system can be described entirely 
by the space coordinates of the vortices and their strengths. Certain invariants of motion 
apply in a two dimensional sense to vortices in free motion, and this leads naturally to the 
concept of a center of vorticity and an effective net vortex strength. 


At some distance from the vortex system the effects of the vortex system become 
indistinguishable from the effects of an equivalent vortex filament of net strength T located 
at the centroid, or center of vorticity, of the system. For the group of discrete vortex 
filaments of strength Tj located behind one half (semi-span) of the wing, the spanwise 
coordinate of the center of vorticity is defined as: 


r ; x: 


X = 


i = 1 


N/2 

Z x i 

i = 1 


(4.15) 


The strength of the equivalent vortex is: 
N/2 

r = Sr, 

i = 1 


(4.16) 


The equivalent vortex filament (EVF) concept was implemented into the VLM wing 
code. The wake is divided in two by a line parallel to the free stream passing through the 
mid-span. Each filament of the equivalent vortex filament pair has a strength equal to the 
net strength of the portion of the truncated wake which it replaces. Thus the net circulation 
is conserved. The filaments of the equivalent vortex pair are placed at the centers of 
vorticity of the divided wake. The wake is divided in two by a line parallel to the free 
stream. Each filament of the equivalent vortex filament pair has a strength equal to the net 
strength of the portion of the truncated wake which it replaces. Thus the net circulation is 
conserved. The filaments of the equivalent vortex pair are placed at the centers of vorticity 
of the divided wake. 

In the standard VLM calculation flow tangency conditions are enforced. 
Aerodynamic performance characteristics are then obtained from an inversion of an 
influence coefficient matrix. This matrix is based entirely on the geometry of the problem. 
The EVF concept was implemented by making this process iterative. The locations and 
strengths of the centroidally placed far-field vortices were iterated upon until flow tangency 
conditions were met. 

Three wake models were considered for the present study. The first model was that 
of the full wake without truncation (namely y t -> infinity), which is the standard way of 
modelling the wake in VLM. The second model truncated the wake at a finite distance 
without replacing the truncated section of the wake. The EVF model also truncates the 
wake at a finite distance, however the lost wake is replaced by an equivalent vortex filament 
pair. 
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Parametric studies were performed to study the validity of the EVF concept. The 
studies were performed on a rectangular wing having aspect ratio of either four or ten, 
with a flat plate airfoil. The planform was at rive degrees angle of attack. The study 
compared overall and sectional load properties. Sectional lift and drag coefficients were 
plotted for various wake truncation lengths. Overall planform lift and drag coefficients are 
plotted as a function of the wake truncation point. Figure 4-5 shows die induced drag 
coefficient of the wing (a quantity that is especially sensitive to the downwash distribution) 
as a function of the truncation distance (expressed in multiples of the chord) for a wing of 
aspect ratio 4. Very small errors occur between the full wake and the EVF model, even for 
truncation at one chord length behind the wing. The most significant errors occur for the 
case with the wake truncated and not replaced with an equivalent far field vortex pair. For 
higher aspect ratio the same trend is observed, but with smaller percentage errors since the 
importance of the wake decreases with increasing aspect ratio. While it might be expected 
that the integrated drag might be relatively insensitive to the wake model, it was surprising 
to rind that die detailed load distribution were also relatively insensitive. Only the case with 
the far wake completely removed showed large deviations in distribution as the truncation 
point was moved closer to the trailing edge. The EVF model showed errors of no more 
than a few percent, and typically less, as the truncation point was brought to within a chord 
of the trailing edge. 

The results reveal a single very important conclusion. The model using an 
equivalent centroidal vortex pair works extremely well, even when truncation of the near 
wake and transition to the pair occurs relatively near the blade. It follows that the rotor 
code should have excellent accuracy in this regard, since it normally uses more than two 
filaments in the free wake, and typically makes the transition from the overlap wake several 
chords behind the blade. 

4.8 Yawed Row 

In forward flight, many sections of the blade are operating with local effective free 
streams that are yawed substantially from the direction normal to die leading edge. Though 
the local flow vector can be substantially affected by wake-induced velocities, the yaw 
angle at any section can be reasonably approximated by 


A = 


= tan'l 


( |icos a s cosy 
+ pcos a s sin\|/ 


(4.15) 


It is clear that for rotors at high advance ratio substantial yaw angles will be present when 
blades are located at azimuth angles 0* or 1 80*, particularly for inboard sections. These 
yaw angles have several important effects on rotor performance. 

First, as noted above, the presence of yaw can substantially increase the allowable 
maximum lift coefficient on a blade section. The lift on a blade section is conventionally 
defined using only the flow component normal to the leading edge as a reference velocity. 


With c^ defined on this basis, then the maximum lift coefficient that can be carried is 
increased to 


( c lmax )a = 0 
cos A 


(4.16) 
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This delays the onset of stall in that larger sectional lift coefficients can be attained before 
the stall model is invoked. In addition, the presence of yawed flow changes the direction 
and magnitude of the profile drag on the section, in that the drag is assumed to act parallel 
to the local free stream and not parallel to the local chord line. 

The presence of yawed flow also has a significant effect on the lift in the vicinity of 
the blade tip. For wings in rectilinear flight or for rotor blades operating at azimuth angle 
90* or 270’, the flow near the tips is ordinarily essentially parallel to the chord line, and the 
wake leaves the trailing edge smoothly. However, in the highly yawed flow that can occur 
in the fourth and first quadrants, the side edge of the blade becomes an effective "trailing 
edge" where a Kutta condition should be applied. Such an effect is not captured in the 
basic vortex lattice model since the trailing legs of the vortex quadrilaterals in the last row 
at the tip overlap one another, so that they effectively constitute a single vortex line at the tip 
(see Fig. 3-9). Once the wing or blade experiences substantial yawed flow, leaving this 
treatment unmodified can lead to errors in the load distribution near the tip. As suggested 
in Figure 4-6 , the quadrilateral edges at the tip of the lattice will ordinarily generate suction 
forces in the presence of yawed flow. In point of fact, in this flow configuration, the side 
edge of the blade should be treated like an extension of the trailing edge, with a condition 
imposed to null the forces carried at this edge. 

This can be accomplished by allowing the lattice trailers to depart the blade parallel 
to the yawed flow instead of parallel to the chordline, as depicted in Figure 4-7 . The tip 
edge trailers are allowed to move with the local flow as the blade moves from azimuth angle 
270 to azim uth angle 90, though they are fixed parallel to the chord over die rest of the 
azimuth. This treatment can be invoked or bypassed as desired, but it typically has the 
effect of increasing the blade loading near the tip relative to the unmodified lattice layout; 
invoking the yawed tip trailers and imposing the Kutta condition on the tip edge eliminates 
the suction force associated with yawed flow and typically causes the entire load 
distribution at the tip to shift slightly and become more positive. 


53 




Figure 4-6. Top view of typical rotor blade 5 x 20 vortex lattice at advance 
ratio 0.4 , azimuth angle 0° . Yawed flow due to 

produces suction (downward force) on side edge in the fourth 
and first quadrants. 



Figure 4-7. Same blade as Figure 4-6 with tip trailers released at local flow 
yaw angle. Kutta condition enforced on blade tip. 
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5. FORMULATION OF THE ROTOR BLADE DYNAMIC MODEL 


The structural properties of the helicopter blade are clearly critical to the evaluation 
of the response and performance of a helicopter in forward flight. Thus, one of the major 
efforts within the development of the RotorCRAFT code has been to incorporate a realistic 
finite element (F.E.) representation of the blade. This section discusses the formulation of 
the dynamic model itself and outlines its capabilities and limitations. The description of its 
coupling with the aerodynamic model discussed above is reserved for Section 6; further 
description of the inputs required for this portion of the analysis are given in Reference 39. 

A preliminary word on the reference frames used in this discussion is in order. 
Figure 5-1 shows the axis systems used in this calculation. The basic reference frame in 
the calculation is the wind axis frame, with X w pointing into the relative wind, 2^ 

pointing down in the plane containing the rotor shaft and Y w at right angles to these axes. 
The shaft axis is reached by a rotation of an angle a s around the Y w axis. The blade 

axis system is then reached by a rotation of an angle (180-\|/) around the Z s axis, so that 
the Z s and Z w are coincident. In the discussion that follows, the blade axis system is 
also referred to as the 'global' axis system to distinguish it from local axis systems attached 
to each of the elements used in the dynamic model. 


5.1 Finite Element Structural Model of the Helicopter Blade 

The particular finite element (F.E.) model used here to represent the helicopter blade 
accounts for extension, twist and transverse bending displacements. To accurately simulate 
these deformations, the blade is discretized into a number of beam finite elements each 
having a total of 14 degrees of freedom (d.o.f.). Stiffness properties for each element are 
computed from the cross-section geometry and material properties supplied by the user. 
Similarly, the blade mass distribution is used to both define elemental mass matrices and 
also to compute the contributions of blade rotation inertia forces to the stiffness matrices 
(geometric stiffening) and nodal forces. The resulting elemental mass and stiffness matrices 
are then assembled and any constrained d.o.f. eliminated to finally yield the corresponding 
global matrices for the complete blade structure. The approach taken is similar to previous 
implementations of F.E. methods for rotorcraft applications, such as Reference 40. 

The resulting mass and stiffness matrices serve as inputs to the standard eigenvalue 
problem which must be solved to obtain the modal frequencies and shapes for the F.E. 
model. It is pointed out that the modal properties are dependent upon the frequency of 
blade rotation since the geometric stiffening is proportional to the square of this frequency. 
The generalized eigenvalue problem is solved by a standard Jacobi iteration technique. 

The transfer of information between the structural and the aerodynamic models in 
RotorCRAFT occurs chiefly via the modal properties of the blade. In essence, the mode 
shapes are used to compute generalized modal forces from the distributed aerodynamic 
forces. These modal forces drive the corresponding modal responses. These responses in 
turn are used in conjunction with the mode shapes to determine the instantaneous 
displacements and velocities at any point along the blade. Finally, the blade deformations 
and deformation rates provide the necessary information to update the flow field and 
aerodynamic forces, as described in the previous section. Hence, one iterates toward a 
steady state solution where the blade displacements and their rates are constant at each 
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a) Rotation of wind axis system about to reach shaft axis system 
(Z s parallel to rotor shaft) 



b) Rotation of shaft axes about Z s to reach blade axis system (X,Y,Z) 
Figure 5- 1 . Axis system definitions. 




Figure 5-2. Correct and incorrect alignment of the elastic axis (E.A.) 
between adjacent elements. 
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azimu thal location. The remainder of this section explains in greater detail the derivation of 
the structural model and of the mode shapes and generalized forces that enter into the 
dynamic computation. 


5.1.1 Assumptions 

The assumptions inherent in the blade model and geometry are stated below: 

• The blade displacements are of sufficiently small magnitude that: 

linear constitutive relations between stress and strain are applicable, 

the transformation matrix relating the local axes of each element may be 
regarded as constant and equal to the corresponding matrix in the 
undeformed state, 

negligible approximation is made in assuming that rotations due to 
deformation commute, and, 

the twist, bending and extension deformations may be linearly 
superimposed. 


• The blade material is assumed isotropic and the stress-strain relation obeys 
Hooke’s law. 

• The elastic axis for each element is defined. The elastic axes of any two adjacent 
elements coincide at their mutual joining section (see Fig. 5-2). In other words, the elastic 
axis is continuous along the blade. This is necessary to correctly define the assembly of the 
individual blade F.E. s. 

• The principal axes of the cross-section for each element are assumed to be 
perpendicular to the elastic axis of that element. This implies that if there are sweep and 
anhedral changes between consecutive elements then their principal axes will not coincide at 
their mutual section. The degree of approximation introduced into the bending calculation 
will increase with the amount of sweep and anhedral change between adjacent elements. It 
will also decrease with slenderness of the element since the discrepancies resulting from 
non-alignment of the principal axes occur locally in the neighborhood of the joining 
section. 

It is pointed out that the last two assumptions are mainly due to the fact that warping 
effects are modelled in the analysis. One of the chief advantages of the finite element 
method is its versatility in the assembly of the constituent elements. For simple elements, 
e.g., pure beam elements and bar elements, one is free to assemble the components in 
whatever orientations one chooses. Furthermore, discontinuities in the mass and stiffness 
properties from element to element are permitted. However, when modeling warping 
deformations the line of shear centers, or elastic axis, plays a significant role. Tfie current 
formulation approximates the elastic axis by a sequence of straight line segments and it is 
the desire to accurately represent the elastic axis that results in the preceding last two 
assumptions. Thus to the extent that warping effects are significant, failure to satisfy the 
last two assumptions and suitably approximate the elastic axis leads to error in the solution. 
In most cases however, and for the closed tubes representative of helicopter rotor blades, 
warping effects will be dominated by deformations arising from pure bending and torsion. 
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and thus violation of these assumptions will not lead to significant error. This has been 
verified by numerical testing of the FE model for loaded structures containing 90° elbow 
joints and discontinuities in die beam stiffness properties. 


5.1.2 Blade Geometry 

Each of the blade segments defined in the RotorCRAFT blade geometry input is 
subdivided into structural finite elements in the manner specified by the user. Tlie global 
axes for the assembled blade are denoted by XYZ corresponding to the blade axes defined 
in Figure 5-1. Local axes, xyz , are defined for each element such that the x-axis 
coincides with the elastic axis, or the line of shear centers, of the element. Axes y and z 
are derived from the transformation applied to the global Y and Z axes as described 
below. The transformation matrix relating the local element axes to the global axes is 
derived from the local segment layout specifications. The segment geometry is specified as 
follows (see Fig. 5-3): 


(1) The planform is first defined. Each blade segment has length, SL , 
along the global X-direction and chord length, c , in the global Y-direction. The sweep, 

A , defines orientation of the quarter chord line for the segment Note that for non-zero 
sweep, the length of the finite element along the quarter chord length differs from the length 
measured along the blade X-axis. If one finite element is associated with each blade 
segment then the element length shall in fact be: 

A = (5-1) 

cos y 

where y is the anhedral (see step 3 below). 


(2) A camber and then a pre-deformation twist gradient are defined over 
each segment This information is not included in the transformation matrix since it is 
judged that effects due to camber and pre-twist upon structural properties can be more 
accurately specified in the information on blade cross-section properties (see Ref. 39). 
Addition of camber would be reflected in the cross-sectional moments of area and pre- 
twisting would affect primarily the orientation of the principal axes. These parameters are 
directly specified in the blade cross-section input file discussed in Reference 39. 


(3) Anhedral is then applied to each segment about an axis parallel to the 
global Y axis and passing through the left hand end (nearest to the rotor hub) of the 
segment. The direction of this rotation is in the negative Y-direction., i.e., positive 

anhedral, y , results in the blade drooping down. 


(4) Finally, collective pitch in the form of a rotation about the global X-axis 
is applied to the assembled structure. 

This sequence of rotations is used to define the transformation matrix relating the 
local axes to the global ones of the RotorCRAFT code. An additional 180° rotation about 
the global X-axis precedes the above rotations since the local finite element z-axis is 
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positive upward whereas the global Z-axis is positive in the downward direction. The 
preceding parameters are supplied in the blade geometry input file. 


From the above sequence of rotations the local and global axes are related by: 



rx 0 
Y ° 
1 Zo 


+ 


T c Ty T A T i80° 



(5.2) 


where Tjgo 0 , T/\ > Ty and T c are the transformation matrices corresponding to the 

180° rotation, sweep, anhedral and collective operations respectively, and the coordinates, 
XqYqZo , are the global coordinates of the origin of the local axes. Here the origin lies on 
the elastic axis at the end of the element nearest the rotor hub. The combined matrix, 


[Trot] ~ T c TyT^ TjgQQ - 


c c s A' s c s t c A " c c c A" s c s 'j®A s c c y 
s c s A +c c s t c A " s c c A +c c s t s A ~Cc£y 


(5.3) 


Note that the rotations due to deformation can also be referred to the global axes using this 
transformation since the deformations are assumed small and the rotations thus commute. 


5.1.3 Element Degrees of Freedom 


The specification of the shape functions and the fourteen degrees of freedom of 
each element is summarized here. Each element has two end nodes and one node at its 
midpoint, as shown in Figure 5-4. The degrees of freedom correspond to translational and 
rotational deformations at these nodes. The deformation of the element at any point is 
estimated by interpolation of the nodal displacements using the shape functions. Let u , v 
and w denote the displacements along the local x , y , and z axes respectively and let 

0 denote the twist deformation about the x axis. Then the generalized displacement 
vector is defined as: 
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(5.4) 


where the subscript (*), x denotes the derivative with respect to the local x axis coordinate 
and A is the element length. Thus qi and q2 refer to the displacement and 
corresponding slope due to bending in the y-direction at the left hand node. The 
corresponding right hand node deformations are q 3 and q 4 , and so forth for the other 
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displacements. Note that the slopes, v, x and w, x , can be regarded as a small positive 
rotation about the local z-axis and a small negative rotation about the local y-axis 
respectively. 

The transverse displacements, v and w , are interpolated using cubic Hermitian 
polynomials as is the common practice in beam finite element formulation. Quadratic 
polynomials are used to interpolate the torsional and axial deformations. This is the 
simplest element interpolation scheme yielding a consistent formulation for coupled torsion- 
bending (Ref. 40). Specifically: 


v(x) = {0 3 ) T j 

, Q4 > 


w(x)= {O 3 } 


r <15 ' 

q 6 
q7 
q8 J 


0(x)= { 02 ) T i qio } 
l qn J 


r q 1 2 1 

u(x)={02} T j qi3 \ 

l qi4 J 


(5.5) 


where, 


{ 0 3 }=< 


1 - 3^ 2 + 2£ 3 
(5-2 5 2 + 5 3 )A 
35 2 - 2 £ 3 
(-5 2 + 5 3 )a 


A 




J 




\ 


- 35 + 2£ 2 ' 
45 - 4 £ 2 > 

-5 + 2$ 2 


(5.6) 


and ^=x/A . The preceding relations may be expressed in compact form as: 


< 


v(x) 'j 
w(x) I 

0(x) 
u(x) J 


[ 0 ]{q} 


(5.7) 


where {q} is the vector of generalized displacements and [ C> ] is a (4 x 14) matrix 
appropriately constructed from Equations 5.5 and 5.6. 

Finally in the formulation of the elemental stiffness and mass matrices it is valuable 
to define principal axes of the cross-section, T) and £ , which are oriented such that: 

JilCdA = 0 (5.8) 


The angle (3 is then the angle between the local y and q axes. 

The global degrees of freedom are obtained by resolving the deformations along the 
global axes using the transformation matrix derived previously. At an end node, all of the 
three translational and three rotational cLo.f. are available (since the slopes of the transverse 
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bending displacements correspond to rotations). Thus, the translation between the local 
element d.o.f.s and the global ones is straightforwardly achieved using the transformation 
matrix, T rot : 



where Rx > R Y » and R Z ♦ are rotations due to deformation about the global XYZ 
axes respectively. At the mid-node the preceding translation is more involved since only 
two d.o.f., the twist and axial deformations, are available in the local axes and additional 
constraints are necessary to uniquely define the twist and stretch in the global directions. 
One approach would be to specify the four remaining local d.o.f at the mid-point by 
inteipolating from the end-nodes using the element shape functions, i.e. evaluating v(A/2) , 
w(A/2) and their slopes using Equations 5.5-5.7. These together with the local twist and 
extension deformations completely define the six local displacements from which the global 
deformations readily follow. However, it was found that this led to numerical problems in 
the resulting transformation matrix, since the complete element is singular for certain blade 
geometries. This might be expected from the observation that 18 global displacements (6 at 
each node) have been defined in terms of only 14 element d.o.f. Hence, the inverse 
transformation from global to local deformations is in fact non-unique. The approach taken 
here is to simply define the global deformations to coincide with the respective local ones at 
the mid-node, i.e., 



(5.10) 


This both simplifies the transformation and results in an orthogonal element transformation 
matrix, i.e., if the elemental transformation which will be composed of elements of [T rot ] 

is denoted by Pgl] so ^ at (qJg = ITgl] (qIl then [Tgl] * = [Tql]^ . An 
alternative procedure would be to eliminate the mid-node d.o.f. using static condensation. 
However, this is unnecessary in light of the small number of d.o.f.s of the fully assembled 
model, the additional programming complexity and the further approximation that would 
thus be introduced. 


5.1.4 Derivation of the Element Strains and Stresses 

In order to compute the elemental stiffness matrix the strains arising from the 
preceding displacements must be evaluated. The nonlinear expressions for the strains are 
stated: 


e xx " 

u, x + ( 'PQjX )>x + 2 ( ^ + C ) ®»x 


+ 

\ v, x 2 - v >xx { T| cos(P+0) - C sin(p+0) ) 


+ 

\ w,x 2 * w >xx ( C cos(|3+0) + T) sin(P+0) } 

(5.11a) 
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Yxrj - 2 e xrj 

('•VC )(e.x+e„i) 

(5.11b) 

YxC = 2e xC 

(^■; + 1)<8.x + ®ol) 

(5.11c) 


and all other strains are assumed zero. Here, 6 n ] is a nonlinear 2 n< * order torsion term, 

and T^x/q.C) is the Saint Venant warping function expressing the out-of-plane 
displacement, u war p , due to torsion: 


u warp( x » T l*C) - '^( x »‘n*C) ®>x (5.12) 

The linear expressions are easily obtained from above: 

e xx = u >x + ( ^6>x )*x 

■ v »xx ( *1 cosP • C sinP } - w, xx { C cosP + T| sinP } (5.13a) 


Yxn - ( C > e.x ; YxC = ( ' + ' n ) e, x (5.i3b,c) 

These strains are expressed in terms of the vector of generalized d.o.f., (q) , and the shape 
functions and their derivatives w.r.t. x , by substituting for the occurrences of u , v , w , 

and 6 using the expressions, Equations 5.11-5.13. This results in: 


-xx 


- [ Tj cosP - Csinp ] {O3 } 
, - [ Ccosp + TjsinP ] {d>3 ) 

¥,x {<*>2) +'¥{* 2 ) 

l {<*2’} 

{ B! ) T (q) 


T 


r {q} 



(5.14a) 


(5.14b,c) 


The corresponding stresses are derived from the Hooke’s Law: 


°xx 


Ee 


xx 


(5.15a) 


Gxi] = G Yxr| ; OxC= g YxC 


(5.15b,c) 
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5.2 Derivation of the Stiffness and Mass Matrices Via Hamilton’s Principle 

The discussion above defines the relationships for stresses and strains for the 
current F.E. formulation. This section deals with the generation of the mass and stiffness 
matrices that permit these elements to be coupled into the equations of motion of the rotor 
blade. There are several ways of obtaining the stiffness and mass matrices. Here the 
Extended Hamilton’s Principle has been applied. One begins by deriving expressions for 
the kinetic energy and the strain energy of each element. Any remaining forces not 
contained in the kinetic and strain energies are accounted for by an external virtual work 
term. Application of the variational calculus to the Lagrangian, £ = (Kinetic Energy) - 
(Strain Energy), yields the equations of motion that govern the blade dynamics. The 
variation of the strain energy in this case is identical to the expression for the internal virtual 
work. Since the internal virtual work can be written down directly, it is unnecessary to 
execute the intermediate step of obtaining an expression for the strain energy. Hence, the 
actual com putati on of the strain energy is circumvented in the analysis to follow. 


5.2.1 Evaluation of the Finite Element Stiffness Matrix 

The structural stiffness matrix K is composed of a sum of three terms Kg , K r , 
and Kj , where Kg is the usual material stiffness matrix obtained from the internal 
virtual work, K r accounts for additional geometric stiffening due to blade rotation, and 
Ky is derived from the variation of the kinetic energy. K r is derived by expressing the 
forces and moments due to rotation explicitly and regarding these as externally applied 
loads that are accounted for by an external virtual work term. This approach is frequently 
used in buckling analysis and bending and torsion problems where axial forces are present. 
The geometrically nonlinear terms in the strain expressions are of importance since axial 
forces are present and result in significant additional stiffening of the blade. If moderate 
deformations are permitted (i.e. moderate rotations but small strains) then the rotational 
forces and moments are sufficiently well specified using only linear geometrical 
considerations. 

The internal virtual work due is given by: 

W 1 = f a xx • 6e xx + o xT1 • 5y xT j + <*x£ * &Yx£ (5.16) 

Substituting for the stresses and breaking up the volume integral: 

wU J ' 

Substituting for the strains using Equation 5.14 and performing the integrations results in: 

W* = {8q} T [K E ]{q} (5.18a) 

where the desired stiffness matrix is 


f Ee xx • 8e xx + Gy xr] • 5y xT1 + Gy x £ • 8y x £ dA } dx (5.17) 
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[K E ] = 


J { Je{Bi}{Bi> 


+ G ( {B 2 }{B 2 } T + {B 3 }{B 3 } T ) dA }dx 

(5.18b) 


The construction of [ K E ] requires a sequence of integrations, the first being an area 
integration over the area of cross-section at a given station, x , along the element, and the 
second being the integral along the length of the element Evaluation of the area integral 
results in expressions containing various properties of cross-section multiplied by the shape 
functions and their derivatives wxt x . These properties include the cross-sectional area, 

moments of area, area centroids relative to the r) and £ axes, and a total of nine integrals 

involving the warping function, S' . The finite element implementation employed in 
RotorCRAFT does not compute these properties, but instead requires that the various 
cross-sectional area integrals be input directly via the blade cross-section input file. The 
finite element code requires that these properties be specified at the end nodes of each 
element and assumes that they vary linearly between the end nodes. The final integration 
along the element length involves products of the shape functions and their derivatives and 
is effected numerically using Gaussian integration. A list of the cross-section area integrals 
required in the analysis is given in Reference 39. 


In order to determine the contributions to the stiffness matrix and the nodal forces 
due to blade rotation one first defines the position vector of a point on the blade in blade 
coordinates. 


R(X,Y,Z) 

= XI+YJ+ZK 


(!) 

f X ° 1 

x ie + u 

= Y 0 + [T rot ]< 

v+T|cos(P+0)-Csin(P+0) . 

l Zq J 

1 w+qsin(P+0)+£cos(J3+0) . 


which in local coordinates is. 


f x ) 

r x o 1 

f x ie + u 

y | = nwf 

Y0 

► + < v+Hcos(P+0)-£sin(|3+0) • 

l z J 

Zo J 

' [ w+T| sin(p+0)+£cos(P+0) . 


(5.19a) 


(5.19b) 


where xj e is the distance along the elastic axis of the element containing the point. Then 
the body force at any point on the blade due to rotation is: 

f = -p(£2K)x {(QK)x R ) (5.20) 

where p is the density of the blade material, and the unit vectors are aligned with the global 
axes. When expressed in the local co-ordinate system of a particular element, this becomes, 
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z 


= ft 2 [T rot ] T ( 


fX 0 

Y 0 

0 



- 1 

0 0 ■ 


x ie + u 

+ 

0 

1 0 

[Trot]" 

v+T)cos(J3+0)-£sin(P+0) ► 


. 0 

0 0 . 


. w+T)sin(P+0)+£cos(|3+0) . 


(5.21) 


The moment about a point on the element elastic axis, Re a = X ea I+ Y ea J +Ze a K , due to 
the rotational force acting on a volume element located at R somewhere on the blade is, 


m = ( R - ) x f dA dx 

The net forces and moments at a point on the elastic axis defined by x^ are: 
R 


(5.22) 


F = 


f(x,y,z) dA ) dX 


(5.23a) 


M = 


Xea 


J { J [R(X,Y,Z) - R ea (X ca ,Y ea ,Z ca )] x f(X,Y,Z) dA } dX 

(5.23b) 


where R is the value of Xg a at the blade tip. The domain of integration extends from 
Xga to the blade tip since the net force and moment vectors vanish at the blade tip. 


The virtual work expression for the net blade rotation forces undergoing virtual 
deformations that accounts for the geometric stiffening effects is stated (see Ref. 40): 

A 

= J {F x (v'8v'+w'5w') + M yy 8(v"0) + M zz 8(w”0) 

+ Q0'8(0') + ^ M xx 8(v"w'-w"v') ) dxi e (5.24) 


Here, 


Q 


X K (y 2 + z 2 ) dA 


(5.25) 


and the terms, F x , M x , My , and M z are simply the local components of the net 
forces and moments due to rotation summed over the portion of blade lying outboard of the 
point xj e on the elastic axis. Note that the first term in Equation 5.24 represents the usual 
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additional stiffening due to an axial force. The virtual work contribution for this first term 
may be viewed as a differential moment arising from structural deformation, 

dM z = F x (v'dx) and dMy = F x (-w'dx) , moving through virtual rotations, 8v' and 

5(-w') . Note also that neglecting the virtual work term. Equation 5.24, will lead to totally 
erroneous results since geometric stiffening plays a major role in the range of angular 
velocities typical for helicopters. In fact, failure to retain of Equation 5.24 results in a 
softening of the blade, which is clearly incorrect. 

Equations 5.19 to 5.24 contain all of the information necessary for the computation 
of [ K r ] . The remaining procedure is laborious, but straightforward and is briefly 
summarized below: 

• Resolve all forces, moments, and position vectors appearing in the above 
equations. Equations 5.19 to 5.24 , in the local element coordinate system. 


• Substitute for all occurrences of u , v , w , and 0 , and their derivatives using 
Equations 5.5 and 5.4. 

• Replace sin(0+0) and cos(J3+0) by the small 0 approximations. 

• Substitute for f , m x , F x , M , and Q in Equation 5.24 and discard all terms 
of order higher than 2. 

• Carry out the cross-section area integrations. As in the computation for [ Ke ] , 
this area integral can be directly expressed in terms of certain cross-section properties. 

Since the blade material density, p , is now present in the analysis these cross-section 
properties will be quantities such as the mass per unit length, cross-sectional center of 
mass, torsional moment of inertia per unit length, etc. The complete list of parameters 
needed is given in Reference 39. 

• Finally, evaluate Equation 5.24. Judicious use of integration by parts where 
possible simplifies the evaluation of Equation 5.25 . For example, the first term. 


A 

x ie 

1 {F x (v'Sv'+w’Sw*) dx ie = 

F x (x ie ) J v'8v'+w'8w'd(J. 


ne - 


x ie=Q 


r 


+ 


<r 


x ie 

f x { Jcv'Sv’+w’Sw’jdp} dx ie 


(5.26) 


The quantity (v'8v'+w'8w') is easily evaluated from Equation 5.7 in terms of the element 
shape functions and the generalized vector of nodal displacements, {q} . Thus the integral 
contained in the brackets {•} can be written down analytically. The final integration along 
the element from xj e =0 to xj e =£ is done by Gaussian integration. 

• The final expression for the virtual work due to blade rotation is of the form: 
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(5.27) 


< = (q) T [K r ]{q) 

2 

where [K r ] is the additional geometric stiffening matrix and is proportional to £2 . 

5.2.2 Derivation of Element Mass Matrix 

The mass matrix is derived from the Kinetic Energy expression for the element. The 
velocity of any point on the blade is a combination of deformation velocities and velocities 
due to die steady blade rotation: 

R = y + flKxR (5.28a) 

where d(*)/9t refers to rates of change occurring relative to axes rotating with the blade. 

■ 

The components of R in the local element coordinates arc: 


Rx 


r ^ 

• 

u 

* Ry 

• 

^ Rz * 

> = < 

v-[T|sin(|5+0)+£cos(P+0)]0 
c R+[t|cos(J3+0)-£sin(P+0)]0 j 


+ [Tro,] T 


0 -£2 0 

r x o i 


x+u 


£2 0 0 

| Y 0 1 

► + [Trot] - 

v+ilcos(f}+0)-£sin(P+0) 

► 

_ 0 0 0 _ 

1 Zo J 

1 

k w+ , nsin(P+0)+^cos(p+0) . 

) 


(5.28b) 


The kinetic energy for the element is, 
A 


KE = 


1 

2 



P ( Rx + Ry + Rz 


} dA 


dx 


(5.29) 


Again the actual computation of the knietic energy is lengthy, but straightforward. One 
replaces the sin(«) and cos(») by their small 0 angle approximations, substitutes for the 
occurrences of u , v , w , and 0 and their time derivatives in Equation 5.28 and inserts 


the resulting expressions for R x , 



and R z into Equation 5 


.29. The time derivatives: 
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r 


< W f = [*]{q} 
6 

^ u J 


(5.30) 


Performing the integration over the cross-sectional area, A , in Equation 5.29 yields an 
expression containing mass properties of the cross-section similar to those cited previously 
in the computation of the stiffness matrix due to blade rotation. The kinetic energy now has 
the following form: 


A 

K.E. = | f (q) T [A2J(q) + n(q) T [All(q) 

+ Q 2 ( { q } T [ A q] { q ) + {q} T {B()} + Q) ) dx (5.31) 

where [A 2 ], [A i] , [Ao] , (Bo) > and Q) are appropriately dimensioned and are 
functions of the element interpolation polynomials and cross-sectional properties. The 
standard application of the variational calculus in Hamilton’s Principle then results in: 

A 

5(K.E.)= - {8q} T [A 2 ]{q} +^([Ai]-[A 1 ] T ){i) -£2 2 ([A 0 ]{q}+{B 0 })dx 

or, SOLE.) = - {Sq} T ( [M] (q) + [G]{q) + [K X ]{q} - (W) ) (5-32) 


A A 


where, [M] = f [A 2 ] dx ; 

[G] = Q J [A!]-[Ai] T dx 

u 

A 

A 

[K T ] = Q 2 f [A 0 ] dx 

and { f rotJ = Q 2 f (Bo) dx 

(5.33) 

u 


[M] is recognized as the mass matrix for the element, [G] is an anti-symmetric gyroscopic 
matrix, [Kj] contributes to the stiffness matrix, and {fjot) represents the nodal forces 
due to blade rotation. It is pointed out that the matrix, [K X ] , and nodal force vector, 
(frot) » could also be obtained by considering the body force due to blade rotation, f , as a 
distributed external force. Imposing virtual displacements upon the blade under this 
distributed load results in the formation of an external work term: 
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(5.34) 


A 



^ f x 8u + fy8v + f z 8w + (yf z - zfy)80 dA dxj e 


Substituting for the forces, fx> fy, h via Equation 5.21, replacing the deformations, u , v , 

w , and 0 by their approximations, Equation 5.5, and evaluating the integral reproduces 
the [Kj] and { f ro t ) derived via the Hamilton’s Principle. 

• 

It is easy to show that the gyroscopic term [G]{q) does no work, i.e., it 
introduces no damping. The elements of [G] are small in comparison to the 
corresponding mass and stiffness matrix components, and their contribution to the overall 
solution is to shift the natural frequencies slightly and introduce a small degree of coupling 
between the modes obtained when neglecting [G] . If the gyroscopic terms were to be 
retained in the eigencalculation for mode shapes and frequencies, this would entail 
substantial increase in the programming complexity involved in the eigencalculation. It was 
judged that the additional computational effort involved in accounting for gyroscopic terms 
is unjustified in light of their relative unimportance in the context of helicopter forward 
flight dynamic analysis. Hence, in RotorCRAFT the gyroscopic terms have been 
neglected. 


5.2.3 Assembly of the Global Mass and Stiffness Matrices 

The element matrices obtained above are employed in the formation of the 
corresponding global matrices for the complete helicopter blade. The assembly process 
involves two sub-procedures: the first involves referring the elemental matrices and nodal 
forces to the global axes, and the second defines the array indexing that relates the local 
degrees of freedom for each element to the global ones. 

Rotation of the element matrices and nodal forces into a global coordinate frame is 
accomplished in the standard manner 


[ Kgjobai ] = [T GL ] [ Kjocjj ] [Tql]T ; [ Mgioba] ] - ITgiJ [ ] [TglJ 

( F™ I global = (F r0t } local [T G lJ T 


as may be easily verified by noting that the kinetic and potential energies and the virtual 
work are independent of the choice of reference frame. Here, PglJ is the transformation 
matrix described in Section 5.2 relating local and global generalized d.o.f.: 
(Xglobal) = [TGLHqiocal) • 

The blade elements are then laid end to end in sequence from blade root to blade tip. 
Global deformations are defined as outlined in Equations 5.9 and 5.10. However, the 
ordering of the global degrees of freedom is different from the local ones, Equation 5.4. 
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Each element degree of freedom is associated with a global one via an indexing array or 
splay matrix, C(k,ie) , where k is the local degree of freedom (k=l,2,...14) , ie is the 
element number, and C(k,ie) is the global degree of freedom. Having specified a C(k,ie) 
for each element then the global matrices may be constructed by ‘splaying’ components of 
the elemental mass and stiffness matrices into their corresponding positions in the global 
matrices. For example, the [i j] entry of the elemental stiffness matrix for finite element, 
ie, is added to the [ C(i,ie),C(j,ie) ] entry of the global stiffness matrix. In like manner, 
the global nodal force vector is built up from nodal forces for each element 

It remains to specify the actual ordering of the global degrees of freedom. Degrees 
of freedom are numbered upwards starting at the blade root Using the definitions for 
global displacements and rotations given in Section 5.2, Equations 5.9 and 5.10, the global 
degrees of freedom are summarized in Table 5-1 (see also Fig. 5-4). 

Thus far, the boundary conditions at the root have not been imposed, and some of 
the root degrees of freedom are eliminated when these are applied. For articulated blades, it 
is implicitly assumed in RotorCRAFT that the blade is freely hinged in both Rap and lag 
directions, but that the remaining degrees of freedom at the root - the three translational 
displacements and the twist about the X axis - are constrained. Thus, the rows and 
columns of the global mass and stiffness matrices corresponding to these four degrees of 
freedom are removed prior to conducting the modal analysis. For cantilevered blades, by 
definition all root deformations are zero and thus all six degrees of freedom at the root are 
similarly eliminated. However, although the rows and columns corresponding to the 
constrained degrees of freedom are removed when computing the natural blade modes, they 
are employed when determining the inertia forces introduced due to cyclic pitching. 

Let the global mass and stiffness matrices for the blade be given by: 



Mcc : Mq- 


Kcc : ^cr 

M = 



; K = 

• • « • • ♦ • « + 


- Mfc • Mjp _ 


. K rc : K n . 


(5.35) 


where the matrices have been partitioned into submatrices associated with constrained and 
unconstrained nodal degrees of freedom respectively. Submatrices, [MJ and [K^] are 
the mass and stiffness matrices for unconstrained deformations and are used in computing 
modal properties. The constrained degrees of freedom are associated with the translation^ 
deformations, U , V , and V , and also the twist rotation, Rx , at the blade root. For 
cantilevered beams, the constrained degrees of freedom also include the rotations, Ry and 
Rz . Of particular interest are the matrix partitions associated with the root twist 

displacement, 0 C , which shall be required for computing the effects of cyclic pitch. 
Although the twist deformation is zero at the root, it is clear that there is a finite rigid body 
rotation occurring at the root due to cyclic pitching and this will introduce an associated 
inertia term into the equations of motion governing the remaining degrees of freedom. This 
is expanded upon in Appendix A. 
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TABLE 5-1 


Specification of the Global Degrees of Freedom For Element, ie. 

Global Deformation Degree of Freedom 


Left-hand node of element, ie: 

V 

08ie-7 


Rz 

08ie-6 


W 

08ie-5 


Ry 

Q8ie-4 


Rx 

08ie-3 


u 

Q8ie-2 

Mid-node of element, ie: 

qio 

08ie-l 


qi3 

08ie 

Right-hand node of element, ie: 

V 

Q8ie+1 

Rz 

08ie+2 


w 

Q8ie+3 


Ry 

08ie+4 


Rx 

08ie+5 


u 

Q8ie+6 
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5.3 Computation of Modal Properties 

The dynamic analysis in the RotorCRAFT code is formulated in terms of the 
response modes of the blade. Thus, the information contained in the blade mass and 
stiffness matrices is used to obtain a set of mode shapes and their corresponding modal 
frequencies and masses. The standard computation of these modes entails solving the 
generalized eigenvalue problem: 

[K rT ]{x i }=0) i 2 [M rr ]{x i } (5.36) 

where coj is the natural frequency of the i-th mode, and {xj} is the corresponding 
eigenvector, or, mode shape. 

Equation 5.36 is solved using the generalized Jacobi iteration method which is 
efficient and accurate, and can be programmed compactly. The scheme, which is described 
in Reference 42, solves for all of the eigenvalues/vectors simultaneously. These are then 
ordered in ascending frequency, and sorted into the dominant type of deformation present 
in each mode: bending (flap and lag), torsion or extension. The number of modes of each 
type retained in the dynamic analysis in RotorCRAFT is specified by the user. 

Appendix A contains a description of the approach used for carrying out the 
computation of the generalized masses and the natural frequencies associated with the rotor 
blade. Additional information on the input required for the RotorCRAFT dynamic model 
can be found in Reference 39. Documentation there and in the Appendix discusses the 
requirements for input to the dynamic analysis as well as the provision for optional 
computation of modal properties given appropriate inputs on the blade's cross-sectional 
properties. 
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6. STRUCTURE AND OPERATION OF THE FORWARD FLIGHT AIRLOADS CODE 


6.1 Outline 

As discussed above, the RotorCRAFT analysis is designed to solve for the 
aerodynamic loads on isolated helicopter rotors in steady forward flight. Presently, the 
rotor hub is assumed to be fixed in space and the rotor shaft is assumed to be oriented at a 
fixed angle with respect to the free stream. Given this and the absence of interference 
effects from bodies or other rotors, the analysis in its present form is most easily 
characterized as a simulation of wind tunnel experiments. This circumstance has 
conditioned the choice of solution method, which is directed at obtaining periodic solutions 
for the wake geometry, the blade motion, and the aerodynamics loads. This objective was 
judged to be appropriate for this particular effort, given the the focus on steady- state 
forward flight. The "transient" solution achieved during the iteration process does not 
represent a time-accurate calculation of rotor aerodynamics; only the converged result can 
be considered physically valid and consistent. The approach used to obtain this converged 
result is now described. 

The outermost loop in the iteration process is the computation of the evolution of 
the wake geometry over one rotor revolution with blade motion and aerodynamic loading 
fixed. The simulation is ordinarily run for a discrete number of rotor revolutions, and 
between each revolution the cyclic pitch inputs and the collective are adjusted to meet the 
conditions specified by the user for thrust and first harmonic flapping response (ordinarily, 
zero first harmonic flapping is desired, though many experimental cases in the literature 
have been run with nonzero flapping). This trim procedure represents the second loop 
within the overall iteration. Finally, the innermost loop is the calculation of the blade 
motion that is consistent with the current values of the control settings; this calculation is 
carried out using the dynamic model described in Section 5 and is performed assuming the 
wake-induced downwash is fixed at the most recent estimate obtained from the outer loop 
involving the wake evolution. 

The first step in starting this overall procedure is to determine an initial estimate for 
the blade motion and the cyclic controls, using a simple strip theory aerodynamic loads 
analysis and an estimate of wake-induced velocity based on a simple, prescribed 
downwash distribution. Using the blade motion and bound circulation estimates so 
generated, the first loop of the analysis proceeds, beginning with a kinematic wake whose 
geometry is determined by the free stream and the downwash at the rotor disk. The wake 
then is allowed to evolve for one rotor revolution, during which the time history of the 
wake-induced flow field at each evaluation point on the rotor blade is computed and stored. 
At die end of this revolution, this time history of wake-induced flow is passed to a trim and 
dynamics routine that updates the control settings and the blade motion solution to be 
consistent with the free wake flow field. 

Using the updated downwash, the program then checks to see if the first harmonic 
components of rigid flapping are within a specified tolerance of the desired level. If not, 
the cyclic pitch is adjusted and the blade motion calculation repeated until they are. The 
iteration history is tracked and previous results are used to accelerate the convergence 
process. Once trim has been achieved, the thrust is checked to see if it lies within a 
specified tolerance of the desired level. If not, the collective is adjusted and the steps just 
above are repeated until this condition is met. This process has proved to be sufficiently 
robust to handle a wide variety of rotor configurations, as will be evident from the sample 
cases presented below. A flow chart of the sequence of events in the trim cycle is given in 


75 



Figure 6-1; more information in the input requirements for the analysis can be found in 
Reference 39. 


6.2 Blade Dynamics and Trim Procedures 

The following discussion supplements the outline just above and provides 
additional background information on the operation of the analysis, in particular methods 
used to obtain the rotor trim. The determination of the initial blade motion and control 
settings (collective and cyclic) to be used with the free wake rotor calculation ordinarily 
proceeds in two steps. First, after input and initialization, the strip theory aerodynamic 
analysis of the rotor blade is used in conjunction with a harmonic balance dynamic analysis 
to determine a preliminary estimate of the blade motion. Hus calculation uses a simple 
inflow downwash model similar to that in from Reference 21 when determining generalized 
forces: 


w(r,\jr) = 



OR (l + VI (£} cosy -pas) 


( 6 . 1 ) 


The harmonic balance solver invokes a blade dynamic model involving only rigid flapping 
and first elastic bending. The calculation is particularly straightforward because of the use 
of approximate mode shapes and because of a simple treatment of compressibility effects 
that allows analytical evaluation of the generalized forces on the blade. Higher bending 
modes and torsional deflection are assumed to be zero for the purpose of this preliminary 
calculation. 

This strip theory/harmonic balance calculation contains an internal trim iteration so 
that the blade motion, bound circulation, and control settings are derived for the specified 
flight condition . These values are used to initialize the appropriate quantities for the full 
free wake/vortex lattice/finite element calculation. It is sometimes desirable to bypass this 
strip theory initialization for repeat calculations at a specified flight condition since the trim 
control settings from a prior run will usually provide a better starting point for the overall 
trim calculation than an initialization based on the simple flow field model described above. 
The implementation of this option is described in Reference 39. 

This strip theory blade dynamics analysis was developed originally for the work 
described in Reference 13 and was retained as the primary source for blade motion 
calculations through the work reported in Reference IS. It has now been superseded by the 
dynamic model described in Section 5, though it may still be invoked to provide blade 
motion information for free wake calculations, if desired. This capability is useful for 
illustrating the effects of the dynamic model of the blade on airload calculations. 

This preliminary solution for blade motion is then used as input to the second 
portion of the initialization procedure. The trim solution derived from strip theory is used 
to drive the vortex lattice model of the rotor blade, but with the wake-induced downwash 
still assumed to be in the form given in Equation 6.1. The use of the vortex lattice model 
will produce new generalized forces, and so a new blade motion and trim solution must be 
obtained to produce a consistent calculation. 
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ISTRM = 2, 


Finite Element Modal Analysis 
to determine the mode shapes, 
frequencies and generalized masses 


ISTRIP = -1 


I stri p = o Strip Theory Analysis 

► to determine an initial 

guess for the blade motion 


Figure 6- 1 . Flow chart of the RotorCRAFT analysis. 
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Figure 6-1 (Corn’d). Flow chart of the RotorCRAFT analysis. 
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To obtain this solution, a call is made to a harmonic analysis solution that is more 
general than the simple harmonic balance solver described above, one that can 
accommodate numerically computed generalized force distributions along the blades. This 
solution method is similar in most respects to that described in Reference 5 . The structural 
displacement of the blade is represented by a summation of orthogonal modes, each 
satisfying the equation of motion: 


d i) 2 *!! 


dr|j 




( 6 . 2 ) 


where Hj is the modal amplitude of the ith mode (r) . The structural displacement is 


z(r,V) = X <t>i(r) Tli(V) i = 2 > N m < 6 - 3 ) 

i 


Here, is the frequency ( normalized by the blade rotation £2 ) of mode <J>j(r) and Gj 
is the corresponding generalized force. The Gj term contains the aerodynamic damping. 
In order to solve the equation using the this approach, damping is required on the left-hand 
side of the equation. An additional artificial damping term , q , is added to both sides of 
the equation, leaving 


9rii 

Fi = Gj + [cj] — ^ 

dy 


(6.4) 


The formulation is simplified by Fourier transforming both the displacements T)i 
and the generalized forces Fj . In the following, the modal subscript 'i' is dropped for 
clarity, but the procedure is applied to each of the modal responses. The solution procedure 
ordinarily starts with an initial guess for the harmonics of the Fourier transformed 

displacements, f| n , based on the strip theory calculation mentioned above. The modal 
amplitude t|j and amplitude rates dqj/dvj/ at the azimuth j are obtained from the inverse 

Fourier transform of the q n 's . Given the values of Tjj and dqj/3\|f , a new value for the 
generalized force Fj is determined at this azimuth using numerical integration of the 
generalized forces over the vortex lattice. The difference between this new value and the 
previous value of Fj creates a difference in each of the harmonics of the Fourier 

transformed generalized forces, AFn , which is used to update the corresponding Tj n 's . 
This process is then repeated at the next azimuth j+1 using the the new T) n 's. After a 
complete revolution, the RMS change in the rj n 's is evaluated over the previous cycle and 
the process is stopped if it is less than 1.0% of the RMS value of the q n 's themselves. 

With the blade dynamics converged, the issue of selecting the pitch controls to 
achieve the desired trim state can be addressed. The options are: 

i) For thrust: fixed collective or fixed Ct with collective adjusting until the 

requested CT is obtained. 
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ii) For first harmonic flapping angles: fixed cyclic inputs Ai w and Bi w or fixed 
first harmonic flapping ais and bis .with Aiw and Bi w adjusting until the selected 
values for ai s and bis are obtained. 

Once the blade motion has been determined, the problem of adjusting cyclic pitch to achieve 
the requested level of flapping is addressed. The partial derivatives 9ais/3Ai w , 
dbis/dAi w , dais/5Biw , and 9bis/5B j w' are detemuned numerically. If the first 
harmonics of the rigid flapping mode differ from the user requested values of ai s and bi s 
by more than a specified tolerance, Aiw and Biw are adjusted according to the 
calculated partial derivatives and the blade dynamics equations are re-solved. This process 
is repeated until the convergence condition has been satisfied. Once this step has been 
completed successfully, the thrust coefficient is checked to see if it differs from the desired 
value by less than a percentage specified by the user. If it does not, the collective is 
adjusted and the steps just described are repeated until it does. The partial derivative 

dCr/900 is calculated while iterating in order to facilitate rapid convergence 

All of the steps described to this point in the section are carried out between 
revolutions of the main rotor wake. The trim process proceeds with a fixed induced 
velocity distribution from the previous cycle of the wake calculation. After the blade 
dynamics calculation is complete, the new blade motion and bound circulation distribution 
around the azimuth are available for application to the next cycle of the free wake 
calculation. The new free wake cycle is performed using updated information on wake- 
induced velocity on the blade and around the azimuth. Then die blade dynamics is repeated 
once again. 

Currently, the user has the option of applying overall convergence criteria to the 
complete calculation, as well. The criteria presently in place allow the calculation to be 
terminated when the thrust coefficient and the cyclic pitch settings are repeatable from one 
free wake revolution to the next within a specified tolerance. Reference 39 contains 
substantial additional discussion of the actual implementation and operation of these 
procedures within the context of the RotorCRAFT code. 
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7. DATA CORRELATION STUDIES 


7.1 Objective 

As is clear from the discussion in Sections 1 and 2, the primary initial motivation 
for the development of the current RotorCRAFT analysis was to capture experimental 
results on the unsteady loading on rotor blades, particularly those that influence vibratory 
loading. In pursuit of this objective, though, it must be kept in mind that the various 
frequency components of rotor loading are closely coupled, and that there is strong 
interaction between the relatively high-frequency components of loading that govern 
vibratory excitation and the low-frequency components that establish blade trim and the 
operating conditions of the rotor. These two problems cannot be easily decoupled and so it 
was judged essential to first address the ability of RotorCRAFT to predict the steady and 
low-frequency loads that establish rotor trim. Indeed, as the extensive discussion by 
Harris in Reference 18 makes clear, even the computation of these quantities is not yet on a 
satisfactory basis for day-to-day calculations in the rotorcraft industry. 

After carrying out correlations focusing on steady and low-frequency loading, 
several correlation exercises pertaining to unsteady loads will be described. These will 
include cases in a wide variety of flight conditions and will examine rotor loading at several 
representative points on the rotor blades in each configuration. These sample calculations 
will provide examples of the wide applicability of the code, but will also note some of the 
limitations inherent in both the analysis and the experimental data . 

In his review of aerodynamic loads data. Hooper (Ref. 6) noted the relative paucity 
of complete and reliable airloads data on modem rotor configurations. The shortcomings 
of existing rotor data sets are also discussed by Johnson in Reference 8. A wide variety of 
problems with existing data are cited by these reviewers, including sparseness of pressure 
transducers for the chordwise integration of loads, incomplete coverage of radial stations, 
and excessive levels of unsteadiness and noise, as well as the presence of poorly 
documented problems with individual experiments. In spite of these reservations, several 
data sets do exist in the public domain that are adequate for illustrating the major features of 
the performance of the RotorCRAFT analysis. Those limitations in the data that may affect 
the conclusions drawn from the correlation will be noted in each case. 

Finally, special attention will be paid to a few particular cases of theory/experiment 
correlation that illustrate the advantages associated with the implementation of full-span 
rotor wake modeling. This will have largely to do with the features of rotor loading that are 
captured with this wake model that are absent from more simplified treatments. 


7.2 Performance and Trim Calculations 

A logical first step in the study of rotor load correlation is the comparison of 
measurements and predictions of integrated rotor loading. A candidate data set for such a 
comparison is the record of wind tunnel tests described Reference 43. These tests involved 
full-scale experiments on both articulated and teetering rotors and measured integrated 
forces and moments for a wide range of shaft angle of attack at advance ratios 0.3 , 0.4 , 
and 0.46 . The comparisons shown here will deal only with the four-bladed, articulated 
rotors within this data set and will cover the following cases and flight conditions: rotors 
with -8* and 0* linear twist, operating at advance ratio 0.3 and 0.4 and at shaft angle of 
attack -5*, 0*, and +5*. The rotor configuration and planform for the -8* case is that of 
the H-34 main rotor, and the 0* twist rotor was apparently identical except for the lack of 
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twist. The planform properties are given in Table 7-1; the structural properties of the blade 
were not specifically stated there but were assumed to be the same as those given for the 
H-34 rotor blade in Reference 44. 

The test results in Reference 43 presented integrated rotor lift and propulsive force 
coefficients normalized by solidity, as well as integrated pitching moment, rolling moment, 
and torque. The shaft axis T , H , and Y forces in the RotorCRAFT output were 
resolved into the wind axis system for comparison to these results. Calculations for a 
range of collective pitch inputs (all trimmed to zero first harmonic flapping, as was 
stipulated in Reference 43) were carried out and the comparison to experimental results are 
shown below. 

The configuration selected for the model was kept relatively simple for these 
calculations. The vortex lattice used on the blades consisted of thirty quadrilaterals 
span wise and a single quad chord wise. The default settings for the core size were selected 
and no adjustment was made during the course of these calculations. The wake structure 
featured a maximum of 16 filaments trailing from Zone 1 and a maximum of six from the 
tip filament region. Zone 2. One turn of full-span free wake was used along with one turn 
of the free wake extensions. Additional far wake was added during numerical experiments, 
but the change in rotor loading due to the addition of rotor wake beyond two turns was 
negligible for the cases examined. NACA 0012 airfoil section data was available in C81 
format for use in providing drag characteristics. The structural and dynamic parameters of 
the blade were input to allow the code to compute the blade's mode shapes and natural 
frequencies. The dynamic model of the blade called for three out-of-plane modes (rigid 
flapping, plus two elastic bending modes). The computation used twenty four time steps 
per revolution and most of the cases discussed below were found to converge well within 
three rotor revolutions. 

Figure 7-1 shows the correlation of rotor lift and torque for the case of advance 
ratio 0.3 , blade twist of -8*, and shaft angles of -5* , 0* , and +5*; the plot shows the 
predicted curve as well as the experimentally measured points. The correlation is close 
over the range of collective pitch settings examined, though the predictions have a tendency 
to overpredict torque over most of this range. The rise in power required that accompanies 
stall appears to be captured, though the power rise begins slightly earlier than the data 
suggests that it should. The fact that the difference between predictions and experiment is 
reasonably consistent across the span for all lift coefficients suggests that profile power 
may be overestimated. One possible source for such an overestimation is the absence of tip 
relief in the drag calculation. At present, there is no allowance for a reduction in the 
effective Mach number in the tip region. Limited numerical experiments with this issue 
have indicated that reasonable approximations to tip effects provide an effect of the right 
sign to account for overprediction of torque though the magnitude may be too small. Also, 
since the largest errors appear around the onset of stall, it may be that the current approach 
to modeling the aerodynamic loads in the presence of blade stall may be inadequate. The 
discussion in the summary below will briefly describe work that is underway to improve 
the stall model. In addition, the sectional drag data from the airfoil tables may be in error. 

Figure 7-2 shows similar results for the case of the rotor with untwisted blades, 
again at advance ratio 0.3. Largely the same results are observed: a close tracking of the 
predicted power coefficient with some small overprediction across the full range of the data 
observed, except for a tendency to underprediction at high thrust levels. Figures 7-3 and 
7-4 show more restricted comparisons of the predicted and measured power over the thrust 
range of interest at advance ratio 0.4 ; the figures show results for the twisted and 
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TABLE 7-1 


SPECIFICATIONS OF THE H-34 MAIN ROTOR 


Number of Blades 

4 

Radius 

28.0 ft 

Chord 

1.366 ft (constant) 

Solidity 

.0622 

Root Cutout 

16% of radius 

Hap Hinge Location 

3.7% of radius 

Twist 

-8* (Linear) 

Airfoil 

NACA 0012 

Rotor Angular Velocity 

Variable in flight tests 
23.2 rad/sec in wind 
tunnel tests 
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Figure 7-3. Predicted and measured performance for a four-bladed rotor 
with -8° linear twist at advance ratio 0.4 and 0° shaft angle 
of attack. 



Figure 7-4. Predicted and measured performance for a four-bladed rotor 
with untwisted blades at advance ratio 0.4 and 0° shaft angle 
of attack. 





untwisted blade, respectively, for the case of cts = 0* . These comparisons show much the 
same trends as the results of the calculations at advance ratio 0.3 . 

Along with integrated performance, the prediction of rotor blade motion is 
important to the understanding of rotor trim. The prediction of the blade s flapping motion 
is in general a more sensitive calculation than the computation of integrated loading. This is 
particularly true at low forward speeds, where the calculation of the effects of the distorted 
rotor wake and its interaction with the blades has a strong influence on the blade motion. 
This sensitivity was well documented in an experimental study using a model rotor (Ref. 
45). This test used a model rotor whose specifications are summarized in Table 7-2. The 
test was carried out in the Boeing Helicopters V/STOL wind tunnel using a four-bladed 
model rotor, the forces and moments on the rotor as well as the blade's coning and flapping 
motion were recorded for a series of runs from advance ratio 0.0 to 0.24 . The runs to be 

discussed here were carried out at essentially a constant thrust level of Cj/o of 0.08 , with 

the shaft angle varied to keep (Xtpp at approximately 1.0*. 

The test was carried out on a balance simulating the rear pylon of a tandem 
helicopter, so the results may well contain body interference effects. The discussion in 
Reference 45 does not elaborate on this possibility, though several data points described as 
being taken in a "previous isolated rotor test " are presented. In addition, the measured 
lateral flapping presented in a nominal hover condition of advance ratio 0.0 was 0.34* , 
suggesting some type of interference effect Finally, the discussion in Reference 45 quotes 
a measurement accuracy of 0.25* on the flapping angles presented. Even though there are 
substantial experimental uncertainties associated with this data, the qualitative results of this 
experiment and their importance to the issues of low speed trim are not in doubt, as is 
discussed at length in Reference 1 1 . 

Several computational efforts (including Refs. 11 and 46) have focused on 
recovering the flapping behavior measured in this test. Previous work in Reference 46 
showed the sensitivity of the predicted flapping to the choice of tip vortex core size with a 
single-tip-vortex wake model. As discussed in Section 3.2.4, one objective of the present 
effort was to reduce the role of vortex core sire as a modelling parameter. The calculations 
presented here for low speed trim all use the default core selection procedure described 
earlier, namely that the core sire of each filament adjusts according to the spacing of the 
filament release points. For these particular calculations, a "floor" or minimum allowable 
core sire of .02R was selected for all filaments. 

For these calculations, a vortex lattice with 30 quadrilaterals spanwise and one 
choidwise was used. The blade dynamics model used a three-mode approximation, with 
the rigid flapping, first elastic bending, and first elastic torsion mode included. The wake 
model featured one turn of full span free wake, three turns of free wake extensions, and 
included the far wake summation described in Section 3. The full-span wake model used a 
maximum of 12 filaments in Zone 1 and two filaments in Zone 2 (the appearance of 
negative tip loading was not anticipated in these calculations and was not, in fact, 
observed). 

Figure 7-5 shows the measured and computed lateral flapping angles for a senes of 
cases from advance ratio 0.06 to 0.24 . Error bars corresponding to the nominal 0.25* 
experimental uncertainty in angle measurement are used with the data points. As is evident, 
the computations fall well within the error bars for the whole range of the test. The peak 
lateral flap amplitude was observed at advance ratio 0.08 , as in the experiment, but the 
trend to reduced flapping at lower speeds could not be reproduced because of poorly 
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TABLE 7-2 


SPECIFICATIONS OF FLAPPING TEST MODEL ROTOR 


Number of Blades 

4 

Radius 

2.73 ft 

Chord 

0.19 ft 

Solidity 

.0892 

Root Cutout 

19.2% of radius 

Flap Hinge Location 

2.3 % of radius 

Twist 

-9.14* (Linear) 

Airfoil 

23010-1.58 

Rotor Angular Velocity 

164.8 rad/sec 
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Figure 7-5. Predicted and measured lateral flap angle vs. advance ratio for 
a four-bladed model rotor at C-p = .007 1 . 


a 1s 
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Figure 7-6. Predicted and measured longitudinal flap angle vs. advance 
ratio for a four-bladed model rotor at Gp = .007 1 . 
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converged time-marching solutions. The lowest speed for which a meaningful answer for 
blade motion could be obtained was advance ratio 0.06 ; this case required nine rotor 
revolutions to converge, while eight were used for advance ratio 0.08 to 0.12 and five 
for all cases with higher advance ratio. The poor convergence properties of low advance 
ratio cases was anticipated and is discussed in Reference 11. With a Lagrangian wake 
model such difficulties are inevitable at low speed, though new methods like that 
documented in Reference 47 may ultimately offer a way to circumvent them. 

Figure 7-6 shows the results for predicted and measured longitudinal flap angle, 
and again the predictions fall within the nominal error bars. The good correlation achieved 
with both lateral and longitudinal flapping at low speed, along with the results presented on 
integrated performance above, indicate that the full-span wake model captures many of the 
important features of the wake-induced flow field for both high and low speed flight. 


7.3 H-34 Flight Test 

Reference 48 contains tabulated measurements of aerodynamic loading on an H-34 
main rotor in flight tests in a wide variety of flight conditions. For purposes of the present 
correlation study, two of these cases were selected. The first was a low forward speed case 
with advance ratio 0.13 and shaft angle of attack -0.7* (Flight 6 of the cases considered in 
Ref. 48). The second case was at a relatively high advance ratio of 0.29 at a shaft angle 
of attack of -6.9* (Flight 19 in Ref. 48). 

The physical description of the rotor blades used here is almost identical to that 
given in Table 7-1. Here, the blades were equipped with a trailing edge tab that was 
deflected 4’ upward between r/R = 0.85 and 0.9 . (The effect of this tab is neglected in 
the calculations that follow). The experimental installation for these tests featured pressure 
transducers placed chordwise along the blade at five to eleven r adial stations from 0.25R to 
0.95R . The time- varying lift force was integrated at each of these stations, and the 
resulting sectional loads are presented in tabular form in Reference 48, along with harmonic 
analysis of these loads and of the pitch and flap motion of the rotor blade. 

The calculations carried out for each of these cases trimmed the rotor using the 
information on gross weight, air density, and tip speed given in Reference 48. The H-34 
used in these tests had a variable RPM rotor, so the tip speed in general changed with each 
flight condition. The tip speed varied for most of the experiments; for the 0.13 case the 
tip speed was 648 fps, while for the 0.29 case it was 690 fps. The thrust coefficient for 
the 0.13 case was 0.0047 , and for the 0.29 case it was 0.0056 . In each of the 
calculations, the rotor was trimmed to zero first harmonic flapping to simulate the zero hub 
moment condition of an aircraft in free flight. Reference 48 notes that the pressure 
transducers on the blade were subject to a 4.5* phase lag; the results below have been 
phase-shifted to eliminate this effect. 

The computations were carried out with a model configuration identical to that 
employed for the full-scale H-34 rotor discussed in the previous section. A 30x1 grid of 
vortex quadrilaterals was used on the blade and a four-mode dynamic model, consisting of 
rigid flapping, the first two elastic bending modes of the blade, and the first elastic torsion 
mode. Both cases used 24 time steps per revolution, though the wake model differed 
substantially for the two cases. For the low speed case, a wake with twelve filaments in 
Zone 1 and two in Zone 2 was used in the full-span wake region, which covered one turn 
of wake. After this point, two turns of two-filament free wake extension were employed, 
followed by two prescribed turns in the intermediate wake and finally the analytical 
summation of the far wake influence. In the high speed case, on the other hand, only two 
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turns of wake were used: one turn of full-span free wake using 16 filaments in Zone 1 and 
six in Zone 2, and a second turn of two-filament free wake extension. It was observed that 
five rotor revolutions were required to converge the solution at advance ratio 0.13 , while 
three sufficed in the 0.29 case because of the convection dominance of the free stream. 
Default core sizes were used in all cases. 

In each case, airloads were computed at four radial stations: r/R = 0.75 , 0.85 , 
0.9 , and 0.95. The results for the 0.13 case are shown in Figure 7-7. Both the data 
and the calculation display the 'peaky' loads in the first and fourth quadrants characteristic 
of wake interactions with the rotor in low speed flight. For this case, the correlation for 
radial stations 0.75 , 0.85 , and 0.9 is quite close, indicating that the wake-induced 
effects in low speed are being captured well. At r/R = 0.95 , however, substantial 
underprediction of loads occurs on the retreating side. The cause of this is unclear, though 
the neglect of the moments caused by the trailing edge tab may be one reason. A trailing 
edge tab deflected upwards would lead to a nose- up pitching moment that could contribute 
to an increase in sectional loading near the tip. 

The results for the advance ratio 0.29 case are shown in Figure 7-8, and they 
exhibit consistently good correlation with the measured loads. However, the loads 
experienced in this particular flight condition have less of the high-frequency loading than 
the wind tunnel tests of an isolated rotor discussed in Reference 44, as will be discussed 
below. 


7.4 SA-349 Flight Test 

The source for the data used in these comparisons is Reference 49, which 
documents flight tests undertaken in 1986 featuring an SA-349. This helicopter has a three 
bladed rotor with blades of radius 17.2 ft. and chord 1.15 ft. Each blade has five degrees 
of washout between the blade cutout and the 92.5% radial station, with an untwisted 
planform outboard of that point. The flight tests covered a wide range of operating 
conditions from advance ratios 0.14 up to 0.38 and included cases in turning flight. The 
level flight case presented here pertains to a rotor at advance ratio 0.14 and a thrust 
coefficient of 0.00427. 

The blades were instrumented with pressure transducers at three radial stations 
(r/R = 0.75 , 0.88 , and 0.97) to provide chordwise load distributions, and these 
pressure distributions were integrated to produce the sectional normal force. Reference 49 
provides filtered frequency spectra of both the normal force and the pitching moment at 
each of the stations, drawn from time histories averaged over six rotor revolutions. As 
discussed by Johnson in Reference 8, the data is subject to considerable unsteadiness and 
noise, and so the repeatability of particular details of the rotor load time history is uncertain. 
In these cases, the primary focus of the correlation study should be on capturing the major 
features of the loading as closely as possible. 

The structural properties of the S A349 rotor blade were taken from the appendices 
of Reference 49, and they were used in the computation of the modal properties of the 
blades. In each of the calculations below, three bending modes (rigid flapping and the first 
two elastic bending modes) as well as the first elastic torsion mode were retained. Default 
core size selection was again invoked. 

The case considered was Test 2 of Reference 49 (advance ratio 0.14, Cj/o = .067). 
These calculations were run with one turn of full-span free wake and two turns of free 
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Figure 7-7. H-34 flight test plot of nondimensional airload vs. azimuth 
advance ratio 0.13 . 
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Figure 7-8 (Cont'd). H-34 flight test plot of nondimensional airload vs. 
azimuth, advance ratio 0.29 . 
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wake extensions; two turns of prescribed wake were also used, and the far field wake 
summation was invoked. The blade loading solutions were found to be adequately 
converged after five rotor revolutions. 

Figure 7-9 shows the predicted and measured time histories of nondimensional 
airload for this case. As is evident, the agreement is good for nearly all azimuthal 
locations, r/R=0.75 and r/R=0.97 . The correlation at r/R=0.88 is worse in certain 
particulars, notably in overpredicting the loads in the vicinity of the 180* azimuthal station. 
It is noteworthy that Johnson in Reference 8 shows a similar overprediction in this region. 
This may suggest that a problem exists with the data, or that there is a shared flaw in the 
analysis of Reference 8 and the current version of RotoiCRAFT. 


7.5 H-34 Wind Tunnel Test 

In Reference 6, Hooper discussed the measured airloads from a wind tunnel test of 
a full-scale H-34 main rotor (Ref. 44). This data set represents one of the few relatively 
complete data sets that are available for correlation studies in the technical literature. As 
Hooper notes, the airload data clearly exhibits some of the dominant mechanisms that lead 
to large vibratory airloads in high-speed forward flight and thus is a very useful point of 
reference for correlation efforts. However, the experimental report indicates that the H-34 
wind tunnel tests were run with a control system that caused the cyclic pitch of two of the 
rotor blades to differ from that applied to the other two. This "split tip-path plane" 
introduces an unwelcome ambiguity into the interpretation of the trim condition for the test 

In addition, examination of the spanwise load distribution in the experimental data 
reveals an unusual feature in the tip loading. Figure 7-10 shows plots of the 
nondimensional thrust measured for the H-34 at the following radial stations: 0.25R , 
0.4R , 0.55R , 0.75R , 0.85R , 0.9R , and 0.95R . Note the presence of a 'notch' in 
the loading distribution at the 0.9R station. The presence of this feature is somewhat 
counterintuitive, since one would typically expect a smoother roll-off in the loading near the 
tip, particularly for azimuth angles around 0* . This feature may be due to a systematic 
experimental problem, such as a biased measurement device or poor calibration. It is also 
possible that it can be attributed to an actual physical phenomenon, such as the irregular tip 
loading due to side edge separation noted in Reference 50. 

The principal physical characteristics of the H-34 blade were given in Table 7-1; the 
rotor tested here featured -8* of linear twist. The rotor's angular velocity for this test was 
23.2 rad/sec , yielding a rotor tip speed of 650 fps , or roughly Mach 0.58 at standard 
sea-level conditions. This tip speed is relatively low by the standards of modem rotorcraft, 
but even so is high enough to produce transonic flow on the advancing tip for flight 
conditions for advance ratios near 0.4 . Appendix B discusses some of the experimental 
uncertainties that can be introduced due to the presence of shock waves in transonic flow 
for such cases. 

Calculations for the H-34 were undertaken using a blade dynamics model that 
included three out-of-plane bending modes (rigid flap mode and two elastic bending modes) 
as well as the first elastic torsion mode. Lag motion was neglected since it was judged 
unlikely that it would contribute significantly to vibratory airloading in this configuration. 
In each case, the computations were trimmed to the measured thrust coefficients and first 
harmonic flapping amplitudes given in the experimental results; in all cases, the nominal 
thrust coefficient was very nearly 0.0037 . For the computations to be discussed in this 
section, a vortex lattice grid using three quadrilaterals choidwise and 30 spanwise was used. 
The default core size option was also invoked, and the calculations were run using the 
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c) r/R = 0.97 
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Figure 7-9 (Corn’d). SA-349 flight test plot of nondimensional airload 

azimuth, advance ratio 0.14 , C T = .00427 . 








option of looking up two-dimensional moment coefficients rather than computing the results 
from the vortex lattice solution itself. 

The computations to be discussed here focus on three cases at advance ratio 0.39 : 
shaft angles of attack -5*, +5* , and 0* . After reviewing the results for airload 
predictions in the first two cases, the results of the 0* case will be examined in detail to 
explore the role of the full-span wake model in the prediction of vibratory airloading. The 
0* case was chosen for particular attention since it involves very close interaction of the 
rotor wake with following blades, and thus the modeling of the rotor wake plays a 
particularly important role in these calculations. 

All of the calculations shown below were performed using 24 time steps per 
revolution, though the predictions were Fourier-decomposed and reconstructed using only 
the first ten harmonics of the rotor rotation frequency. This was done to match the 
harmonic resolution of the experimental data. Also, the measurements include a five degree 
phase lag attributable to the pressure measurement system; the plots of the measured data 
that follow have not b een adjusted to allow for this lag. 

Figure 7-11 shows the measured and predicted time histories of the nondimensional 
thrust loading at five radial stations for the -5* case. At r/R = 0.4 and 0.55 , the 
agreement is very close, with the only significant differences appearing in the vicinity of 

V = 0* . The results at r/R = 0.75 show a significant phase error, though the levels are 
well predicted; part of the phase error can be accounted for from the experimental lag just 
described, but it does appear that the prediction leads the measurement nonetheless. At 
r/R = 0.85 , the level is significantly off and a phase error persists, though the level close to 
the tip ( r/R =0.95 ) is well captured. 

Largely similar results appear in Figure 7-12, depicting the correlation at +5 shaft 
angle of attack. The results at the inboard sections shown in 7-12a and 7-12b are well 
predicted, though more significant errors appear at r/R = 0.75 . Here, the size of an "up- 

down" pulse in the vicinity of \|/ = 0* is overpredicted and the phasing of the loading on 
the advancing side is in considerable error. However, both features improve for stations 

nearer the tip, as shown in Figures 7-1 2d and 7-12e; the size of the pulse in the \y = 0' 
region is closer to experiment and the phasing of the loading becomes considerably closer. 
Figure 7-13 shows a similar pattern for the loading in the case of 0* shaft angle of attack; 
close correlation inboard and near the tip, with more significant errors in level and phasing 
for intermediate stations. 


7.6 Contributions of Full-Span Modeling in Airload Calculations 

These results indicated that many of the features of the airloading are being captured 
for the H-34 case. A topic of special interest at this point is the source of the aerodynamic 
loading that contributes to vibratory loads on the rotor. This section focuses on one 
particular case and discusses in detail the influence of the current full-span wake model and 
the vortex lattice blade model on the prediction of the components of the rotor loading that 
contribute to vibration. The H-34 is a four-bladed rotor, and so the 3P , 4P , and 5P 
(three-, four-, and five-per-revolution) components of the aerodynamic loading will be the 
primary contributors to vibratory loads at the rotor hub. The first priority in this 
discussion is to identify the physical mechanisms leading to such loading. 
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b) r/R = .55 Azimuth Angle (deg) 

Figure 7-11. Nondimensional sectional thrust vs. azimuth angle for the 
H-34 rotor: advance ratio 0.39 , -5° shaft angle of attack. 






c) r/R = .75 Azimuth Angle (deg) 
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d) r/R = .85 Azimuth Angle (deg) 

Figure 7-11 (Cont'd). Nondimensional sectional thrust vs. azimuth angle 

for the H-34 rotor: advance ratio 0.39 , -5° shaft 
angle of attack. 
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Figure 7-11 (Cont'd). Nondimensional sectional thrust vs. azimuth angle 

for the H-34 rotor: advance ratio 0.39 , -5° shaft 
angle of attack. 
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b) r/R = .55 Azimuth Angle (deg) 

Figure 7-12. Nondimensional sectional thrust vs. azimuth angle for the 
H-34 rotor: advance ratio 0.39 , +5° shaft angle of attack. 
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d) r/R = .85 Azimuth Angle (deg) 

Figure 7-12 (Cont’d). Nondimensional sectional thrust vs. azimuth angle 

for the H- 34 rotor: advance ratio 0.39 , +5° shaft 
angle of attack. 





Figure 7-12 (Cont'd). Nondimensional sectional thrust vs. azimuth angle 

for the H-34 rotor: advance ratio 0.39 , +5° shaft 
angle of attack. 
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d) r/R = 

Figure 7-13 (Cont'd). 






Figure 7-13 (Cont'd). Nondimensional sectional thrust vs. azimuth angle 

for the H-34 rotor: advance ratio 0.39 , 0° shaft 
angle of attack. 
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First note Figure 7-14, which shows the top view of the wake structure for the 
H-34 rotor and points out the major features of the vortex wake, including several closed 
loops corresponding to the wake of the negatively loaded tip vortex zone as well as the 
inboard wake, composed of several filaments trailing from the inboard portion of the 
blade. While providing substantial details of the structure of the wake, this figure only 
shows the wake of one blade at one instant and so is of limited usefulness in understanding 
the time varying flow field that characterizes blade/wake encounters. Figure 7-15 does not 
contain the complete details of the full-span wake, but it does show the major features of 
the wake of all of the blades as they interact with die reference blade (shaded black). The 
top portion of Figure 7-15 shows the time history of the upwash velocity (positive opposite 
the shaft Z axis) at the r/R=0.95 station on the reference blade for the shaft angle O' case 
described in Figure 7-13. The letters A-D correspond to the various major wake/blade 
encounters around the azimuth. 

Event A corresponds to the passage of the tip of the reference blade over the tip 
vortex trailed from a previous blade when passing through the third quadrant Event B 
occurs because of the strong upwash/downwash experienced by the blade as it passes over 
the root vortex system. At the time indicated by Event C, the downwash induced by the 
combined effect of the root vortex and tip vortex systems at the reference blade is beginning 
to diminish, gradually reversing into an upwash velocity as the blade tip encounters the 
flow field generated by the opposite sign vorticity that trails from the negatively loaded tip 
region (Event D); note the the existence of negative loading at the tip in Figure 7-13e. 

Because no flow field measurements are available for correlation to the predictions 
given in Figure 7-15, their validity can only be inferred through comparison to the loading 
data. Comparing the time history of upwash to the complete loading time history in Figure 
7-13e suggests the wake-induced loading does behave much as indicated in Figure 7-15. 
However, the relationships involved are still clearer if only the higher harmonics of the 
airloads are plotted. Figure 7-16 shows both the predicted and measured time histories of 
airload at r/R = 0.95 with the lower harmonics (0,1, and 2 ) removed. These residual 
higher harmonics are those that would generate the vibratory loads transmitted to the 
fuselage. Despite the complication introduced by the participation of structural deflection in 
the higher harmonic loading, the correspondence of the loading events to the wake-induced 
upwash is clear. 

The predicted loading associated with event A in Figure 7-15 shows up in Figure 
7-16, but it appears that the analysis overpredicts the size of this event. The characteristic 
up/down loading on the advancing side also is present, though its phasing leads the data by 
somewhat more than the five degrees of phase lag that can be accounted for in the data. 
The level of loading in this event is well predicted, an encouraging circumstance since 
many of the analyses discussed by Hooper in Reference 6 failed to achieve even qualitative 
correlation with the higher harmonic loading in this data set 

Another important feature of the higher harmonic loading is captured, namely the 
large up-down pulse around azimuth angle 0. Figure 7-15 makes clear that this event is 
associated with the passage of the blade over the portion of the wake composed of inboard 
trailing filaments. The presence of this detailed model of the inboard wake is thus clearly 
important for accurate reconstruction of the higher harmonic loading in this case. The 
importance of this mechanism suggests that the inboard wake deserves more detailed 
attention than it has received in rotor wake analyses to date. Analyses that use a single free 
tip vortex and smear out the inboard wake into a few large-core vortex filaments will miss 
this feature entirely, and it is clear that this interaction event makes a significant contribution 
to higher harmonic loading. 
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Closed loops produced by negative 
tip loading on the advancing side 
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Concentrated root vortex from inboard bound 
circulation in the first and second quadrants 


Figure 7-14. Top view of the wake of one blade of the H-34 rotor, 
reference blade at azimuth angle 180° . One and one quarter 
turns of full-span free wake shown. 
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A) Downwash from previous tip vortex 



B) Crossing root vortices 



C) Downwash from root and tip vortices D) Upwash from the tip vortex system 
Figure 7-15. Plot of upwash as a function of azimuth for the H-34 at advance ratio 
0.39 , 0 shaft angle of attack. Events A-D correspond to the indicated 


blade/wake interaction events at \j/ = 315° , 0° , 45° , and 90° 


112 








Figure 7-16. Higher harmonics (3P-10P) of airloading for the H-34 main 
rotor, advance ratio 0.39 , 0° shaft angle of attack. 
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In order to understand the mechanisms involved in these wake blade interactions 
more clearly, additional calculations were carried out with the intention of highlighting the 
contributions of various parts of the wake to the higher harmonic loading. This was 
accomplished by artificially inflating the vortex core size of filaments trailing from various 
portions of the blade. The results of this exercise are summarized in Figure 7-17. In 
Figure 7- 17a, the baseline calculation of higher harmonics airloads for the 0* shaft angle 
case are shown along with the results of two supplementary calculations that were carried 
out with the core sizes for the trailing vortices set to one rotor radius, which should 
effectively delete their role in the calculation. The first case (labeled "tip vortices deleted") 
involved using this core adjustment on the vortices trailing from Zone 2 as well as the 
filaments trailing from the outer half of Zone 1, while leaving the remaining cores at the 
default sizes; the second case ("root vortices deleted") uses the adjusted cores on the 
filament on the inner half of Zone 1, while using default values for the remainder. In point 
of fact, since filaments trailed from the tip zone often loop around to join up with filaments 
from the root zone, this is not a truly consistent procedure for fully deleting the 
contributions of portions of the wake, but it does serve to emphasize the role of particular 
regions. 

Figure 7- 17a shows the results of these core modifications; 7- 17b repeats them 
without the baseline calculation being plotted for clarity of comparison. The deletion of 
root vortex effects removes most of the up-down pulse around azimuth angle 0, while 

leaving strong contributions from the loading events around y = 90 . This is intuitively 
reasonable, since Figure 7-15 indicates that the latter events are driven largely by 
interactions of the blade with wake trailed from the outboard regions of the blade. 

Conversely, the deletion of tip vortex events leaves the interaction around tjf = 0 as the 
major loading event; some residual loading in the first quadrant also remains. This 
comparison indicates the distinctive role of the inboard wake in the full-span wake model, 
and also highlights the contributions of the wake of the tip region. Clearly, both of these 
regions can make important contributions to the higher harmonic airloading of the rotor in 
high-speed flight conditions. 


114 





b) Case with tip vortices deleted and case with 
root vortices deleted 

Figure 7-17. Comparison of predicted higher harmonic airloads using the complete wake 
model with two modified models: tip vortices deleted and root vortices 

deleted. (H-34 at advance ratio 0.39 , a s = 0° , r/R = 0.95 ) 
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8. SUMMARY 


One of the primary objectives of this report has been to motivate the development of 
a new approach to helicopter rotor free wake modeling and to describe its implementation in 
the RotorCRAFT analysis. This description has included a discussion of the principles of 
full-span free wake modeling, focusing on the current approach to the discretization of the 
wake of the rotor blade into a set of constant-strength vortex filaments whose spacing and 
structure visibly reflect the time history of the loading on the rotor blades. As has been 
shown, this approach constitutes an exceptionally general and flexible model for the rotor 
wake of helicopters in both high-speed and low- speed flight 

Along with the description of the fundamental full-span model, another topic of 
interest has been the extension of this basic approach to include the potentially important 
effects of the far wake without impairing computational efficiency. 'Die approach taken to 
this issue has featured the combination of a freely distorting two-filament extension of the 
full-span wake, several prescribed turns of discrete filamentary wake, and an analytical 
summation of remaining far wake effects for completeness. This combination of models 
will allow users of the RotorCRAFT analysis to tailor the structure of the far wake model to 
fit the particular flight condition of interest and to adapt to constraints on CPU time. 
Another important innovation in this last regard is the implementation of vortex elements 
based on Analytical Numerical Matching (ANM). Though not yet fully tested, sample 
problems run to date have indicated that flee wake calculation using ANM may be expected 
to run in 30 to 50 percent less time than calculations using BCVEs with little or no sacrifice 
in accuracy. The ANM wake model may be invoked in the current code as an option, 
though the use of BCVEs is the preferred approach for the present 

The inclusion of a vortex lattice model of the blade has also been described, as have 
the modifications undertaken to the basic lattice model to incorporate the effects of viscous 
drag, yawed flow, and stall. The stall model currently in place has been found to produce 
qualitatively reasonable behavior in regions of high lift and is judged to be a substantial 
improvement relative to unmodified lattice models, which allow arbitrarily high generation 
of lift As noted in Section 4, the current stall model directly replaces segments of the 
vortex lattice in stall with segments loaded according to two-dimensional tabulated 
coefficients. This approach can lead to sharp gradients in load at the boundary of stalled 
regions and so can produce unrealistic load distributions in some cases. An improved stall 
model is currently being explored that will operate by directly modifying the blade-on-blade 
influence coefficients in the vortex lattice solution method to reproduce the nonlinear lift 
behavior of particular sections. This improved model should produce smoother and hence 
more credible load distributions at stall boundaries and, if successful, will be included in 
the follow-on version of RotorCRAFT that is expected in early 1991. 

A substantial portion of the discussion here has involved a description of the finite 
element structural model of the blade. As is clear from the discussion, this implementation 
permits considerable flexibility in the layout of the blade and introduces an important 
increment in sophistication of the model of blade deformation that was absent from earlier 
versions of the RotorCRAFT analysis. 

The final portion of this report, containing the discussion of correlation runs with 
experimental data on rotor loading, was intended to both provide evidence of the code's 
suitability for application to challenging problems in rotor load prediction and to identify 
areas for further development. The calculations of integrated loading indicate a tendency to 
overpredict rotor torque required for a specified rotor lift, and it is speculated that one 
source of this overprediction is the absence of a compressibility tip relief model. This 
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effect and other candidate improvements to the current performance analysis will be 
considered as part of a continuing effort to improve the existing code. 

As to the predictions of unsteady loading, the current version of RotorCRAFT 
produced good correlation with the SA-349 data and the H-34 flight test data, with some 
exceptions at particular radial stations. This indicates that RotorCRAFT performs well for 
rotors at low to moderate advance ratios. Correlation to high speed runs using the H-34 
wind tunnel data set also produced generally favorable results both for overall loading at 
inboard stations and for stations near the rotor tip. However, substantial errors in level and 
phasing were observed, particularly for the .75R and .85R radial stations. The source of 
these errors is a topic of continuing investigation, though the exact extent of the problem is 
difficult to judge, given the uncertainties in the airload data measurements discussed in both 
Section 7.5 and Appendix B. 

The final discussion of higher harmonic loading makes it clear that RotorCRAFT 
can capture significant features of the vibratory airloading of the H-34, and that the use of 
the full-span wake model is an essential component of this capability. This confirms the 
fundamental premise of this work, namely that the use of a full-span wake model is an 
essential prerequisite for accurate reconstruction of vibratory airloads in forward flight. 
However, it is equally clear that other features of the environment of typical rotors, notably 
the presence of transonic flow on the advancing side, can also play an important role in the 
determination of unsteady loading at high forward speeds. It is anticipated that this topic, 
along with enhancement of the current features of the code described above, will be 
prominent among the topics addressed in the continuing development of RotorCRAFT. 
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APPENDIX A 


EIGENCALCULATION SOLUTION METHOD FOR ROTOR BLADE MODAL 

PROPERTIES 


A.l Formulation 

The eigenvectors, {xj} , specify the nodal displacements that characterize a given 
mode shape. It is convenient to interpolate these nodal displacements using the individual 
element shape functions to obtain: 


Oi(r,k) 

where, i 
r 

k=l,2,3 

k=4,5,6 


refers to the particular mode 

radial distance along X-axis of a point cm the blade 

displacements, U, V, -W in the X,Y,and -Z directions 

rotations, Rx, Ry and Rz, about the X,Y,Z axes, or equivalently, 

twist, 0x, about X, and slopes, -3W/9X and 9V/9X, respectively 


XYZ are the blade-fixed coordinates as shown in Figure 5-1. Note the negative signs on 
the W displacement. Positive W occurs in the positive Z direction. However Oj 
should contain -W as the third coordinate. The user of the RotorCRAFT code has the 
option of specifying the mode shapes explicitly rather than computing these from the blade 

stiffness and mass properties. This amounts to evaluating the mode shapes, d>j(r,k), at a 
sequence of radial locations, r=r(ir) , to define an array: 


SHAPES (i,ir,k) = 0>i( r(ir),k ) (A.l) 

The generalized masses are defined as: 

GM(i) = {xi^tMn-] { Xi } (A.2a) 

(in many eigencalculations GM(i) = 1 by virtue of the way that the eigenvectors are 
normalized). If assumed modeshapes have been used then the generalized mass can also be 
expressed as: 


R 

GM(i) = Jc G (r) 


{ <D?(r,l) + <D?(r,2) 44>?(r,3) } + l£(r) 0?(r,4) dr 


(A.2b) 


where 0 ( 3 (r) is the mass per unit length along X and JqW is the second moment of 
inertia per unit length about an axis parallel to the blade X-axis and intersecting the local 
elastic axis. (Note that oq and Jq are referred to the global axes, XYZ , whereas the 

cross-section input parameters in the finite element (F.E.) analysis, o and J m are referred 
to the local axes.) 
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A.2 Rotation Due to Cyclic Pitch 


The cyclic pitch acceleration introduces a forcing term into the modal equations 
involving torsional deformations. In order to account for the cyclic pitch at the root 
additional information must be furnished. 

Let the vector, (r 0 c ) be the vector describing the rigid body displacement of the 

structure due to a unit twist at the root. This can be either computed directly from the 
geometry of the blade, or from: 

{r 0 c } = - [Kjt ]' 1 [Kre c ] (A.3) 

where [K r Q c ] is the column of the stiffness matrix associated with the twist degree of 
freedom of the root node. This is contained in the matrix partition, [K rc ] , defined in 
Equation 5.35. Then the displacement at each node due to rigid body motion arising from 
cyclic pitch is simply, {r 0 c } 0 c (t) . Note that the components of the vector, {r 0 c } , can 
also be obtained directly from the geometry as: 

U = 0 ; V = -i 2 ; W = ry 

®x = 1 ; Ry=o ; Rz=o (A.4) 

where ry and rz are the perpendicular distances of the elastic axis from the cyclic pitch 
axis, X , in the Y and Z directions. 


A.3 Modal Equations 

The remaining derivation assumes a F.E. setting. However, those quantities 
required as inputs to the code for users desiring to input their own mode shapes are defined 
independently of F.E. parameters, as is seen below. The equations of motion defining the 
evolution of the DOF, q^t) , are then, 


[ Mj-0 c : Mjt ] 


r 

• ■ 



® C 


e c 

1 

> + [ Kr0 c : Kir 1 * 

► 

m m 

>• Q x j 


. <lr . 


{F aer°} + {p rotj 


(A.5) 


where [Mj 0 c ] and [Kj 0 c ] are the column of the mass and stiffness matrices associated 

with the root twist, {F aero } are the nodal forces due to aerodynamic loads and { F rot } 
the nodal forces associated with the blade rotation. Note that none of the other constrained 
degrees of freedom, U , V , W , etc., at the root enter Equation A.5 since there is no 

forcing along these directions. The only non-zero root displacements is 0 C . Furthermore, 
it is clear that the generalized coordinate vector, {qj.} , includes both deformation effects 
and the rigid body rotation due to cyclic pitch. Rearranging Equation A.5: 
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[Md ( q r ) + tKn-l (%) = (F aer0 ) + (F rot ) - ( [Mre<J ' e 'c + [Kre c l 9 C ) (A.6) 

The elastic deformation variables {q} are determined by subtracting out the rigid body 
motion arising from cyclic pitch: 

{q} = {qr} - {r0 c }0 c ; (*4) = {<ir> * ( r 6 c }®c ( A - 7 > 

Then, Equation A.6 becomes. 


[MrrKqJ+tK^Kq} = {F aero } + {F rot } 

- ( [M^MMrrHre^ )0 C + ( [K rec ]+[K rr ](r 0c }) 0 C 
= {F aero } + {F rot } - ( [Mj0 c ] + [M rr ]{r 0c } )0 C (A.8) 

where expression Equation A.3 has been used 


Modal coordinates, {rj } , defined by {q}=[X]{t|} , are now introduced. Here the 
modal matrix, [X]=[{ xj}] , is composed of the first N m eigenvectors, {xj} , that have 
been retained in the analysis. These eigenvectors and the corresponding natural 
frequencies, coj^ are the non-trivial solutions to the generalized eigenvalue problem. 
Equation 5.36. Eauation A.8 may then be transformed into the modal variables by pre- 
multiplying by [X] 1 and substituting {q}=[X]{ri } : 

[X] T [ Mj-j- ][X] {il } + [X] T [ Kn- ] [X] { T| } 

= [X] T [ { F aero } + {F rot } - ( [M r0( .] + [Mn.]{r 0c } ) 0 C ] 
or, GM(i) [ TJi + coi 2 r\[ ] = F AER0 (i) + F R0T (i) - MRC(i) 0 C (A.9) 


The generalized modal force due to aerodynamic loading, 
r=RADIUS 


F AERO (i) = 


Hlli 


^ AERO AERO AERO AERO, 

{fy fy fy m x } 


GE 


<I>i(r,l) 

<I>i(r,2) 

-0»i(r,3) 

Oi(r,4) 


r dr 


(A. 10a) 


where f , and f^^ R ^ are the global axes components of the 

aerodynamic force per unit length, and na^p^^ is the aerodynamic moment per unit 

length about an axis parallel to the X-axis and passing through the local elastic axis 
pAERO(i) ? j s approximated in RotorCRAFT via the discrete sum. 
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AFRO f SHAPES(i,ir,2) ] 

F AHKU(i) = y {DDIS(ir) TDIS(ir) MDIS(ir)} { SHAPES(i,ir,3) \ (A.IOb) 
rr l SHAPES(i,ir,4) J 


where, MDIS(ir) = QMDIS(ir) - TDIS(ir)*ELOFF(ir) and DDIS(ir) , TDIS(ir) , and 
QMDIS(ir) , are the arrays used in RotorCRAFT to store the drag, thrust, and pitching 
moment of a chord segment at point, ir (see Ref.39). ELOFF(ir) is the array used in 
RotorCRAFT specifying the distance of the elastic axis behind the quarter chord point in 
the global Y-direction at radial station ir . The modal force due to blade rotation, 

F R ° T (i) = {xi) T {F rotj (A.lla) 


If mode shapes are defined by the user, then 

FROT (i) . ***** {f ROT jROT ,ROT 

r=HINGE 


ROT 

X 


r <I>i(r,l) 

< <&i(r,2) 

-<l>i(r,3) 
. Oj(r,4) 


* dr 
(A. lib) 


where f^^ , fy^T , and are the global axes components of the force per unit 

ROT 

length due to blade rotation, and m^ is the blade rotation twisting moment per unit 
length about an axis parallel to X and coincident with the local elastic axis. Note here the 
sign of d>j(r,3) which is due to the fact that the displacement, d>i(r,3)=-W , is in the 
negative Z-direction. The quantity incorporating the inertial effects of cyclic pitch into the 
modal equation, 

MRC(i) = {xi} T ( [M^ c ] + [M rr ]{r 6c ) ) (A.12a) 

or, if mode shapes are input separately, 


r=RADIUS 

MRC(i) = r { - OQ(r)rz OG(r)rv 
r=HINGE 




<l>i(r,2) 

-Oi(r.3) 

Oi(r,4) 


r dr 


(A. 12b) 


where oq and Jq have already been defined for Equation A.2b and the moment arms, 
ry and i£ , are defined in Equations A.4. 


Finally by dividing through by GM(i) and normalizing the time by x=t/Q , one 
obtains the modal equations: 

Tli" + tO n 2 tl i = — j- q2 [F AER0 (i) + p ROT (i) - MRC(i) 0 C ] (A.13) 

where the normalized angular frequency, o>q = coj/Q . 
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APPENDIX B 


ACCURACY OF MEASURED LOADING ON BLADE 
SECTIONS IN TRANSONIC FLOW 

Assessment of the overall predictive capability of the RotorCRAFT code has been 
primarily based on comparisons of measured airloads with code predictions. These 
airloads were obtained by measuring the pressures at discrete chordwise locations for a 
specified radial station with pressure transducers and then integrating over the chord. As a 
result of spatial constraints, these chordwise pressure transducer locations are often limited 
to relatively small numbers (between five and eleven in a typical test, such as Ref. 44) and 
it is the purpose of this limited study to estimate the error in airloads associated with 
integrating these discrete pressures over the blade chord. 

Transonic airfoil pressure distributions have been measured extensively on NACA 
0012 airfoils in several carefully controlled tests (Ref. B-l). Characteristic of transonic 
pressure distributions on this airfoil is a strong shock wave which leads to pressure 
recovery on the airfoil as the flow approaches the trailing edge. Data is reproduced in 
Figure B-l which gives the measured pressure distribution at Mach .747 angle of attack of 
about 3*. This Mach number and angle of attack were judged to be reasonably 
representative of the blade operating condition at stations near the advancing tip of rotors in 
high speed flight. Pressure taps used to obtain this data in Reference B-l were on the 
average, at a spacing of 3% of chord, much higher than is typically used in helicopter flight 
or wind tunnel tests. 

Since substantial use was made of the Sikorsky H-34 full-scale data, an estimate of 
the error in computed airload was undertaken in an attempt to quantify what constitutes 
acceptable data / rotorcraft airloads comparisons. Table B-l below shows the chordwise 
locations of pressure taps as a function of radial position on the blade. It is significant that 
the tip of the blade is sparsely instrumented and that the highest chordwise resolution of 
pressure is located at 85% of rotor radius. 

Figure B-2 shows the results of integrating the pressure data of Figure B-l as a 
function of discretization of the load distribution. The error presented here is the deviation 
of total sectional load computed with various discretizations of the chordwise load 
distribution to the total load measured by the wind tunnel balance in Reference B-l. The 
error denoted by plus signs results when the load is uniformly discretized in the chord 
direction, while the triangular symbols denote the discretization resulting from the pressure 
tap locations in Table B-l. Note that integration of the chordwise pressure distribution can 
introduce significant errors in the sectional load. These errors will in general be greater still 
for largo* angle of attack or for higher Mach number. It is reasonably safe to conclude that 
for higher Mach number data outboard of the 85% radial position, errors in excess of ± 
8% may exist in the integrated data. 

Reference: 

B-l. Thibert, JJ., Grandjacques, M. and Ohman, L.H.: "NACA 0012 Airfoil" AGARD- 
AR-138, Report of the Fluid Dynamics Panel Working Group 4, May 1979. 
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Figure B-l. Chordwise pressure distribution over a NACA 0012 airfoil. 




TABLE B-l 


Location of Pressure Taps (Ref. 44) 
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Figure B-2. Percent variation of section load (from integration) divided by 
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